ViewVC Help

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

root/group/branches/no-template-branch/OOPSE/libmdtools/ZConstraint.cpp

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 671 by mmeineke, Fri Aug 8 17:48:44 2003 UTC vs.
Revision 736 by tim, Thu Aug 28 21:09:47 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	<
55	<	indexOfAllZConsMols = index->getIndexData();
56	<
57	<	//the maximum value of index is the last one(we sorted the index data in SimSetup.cpp)
58	<	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){
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 error: index is out of range\n");
61	<	painCave.isFatal = 1;
62	<	simError();
63	<
64	<	}
65	<
60	>	"ZConstraint Warning: unknown force substraction policy, "
61	>	"PolicyByMass is used\n");
62	>	painCave.isFatal = 0;
63	>	simError();
64	>	forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
65	>	}
66		}
60	–
67		}
68
63	–	//retrive sample time of z-contraint
64	–	data = info->getProperty("zconstime");
69
70	+	//retrieve sample time of z-contraint
71	+	data = info->getProperty(ZCONSTIME_ID);
72	+
73		if(!data) {
74
75		sprintf( painCave.errMsg,
#	Line 89 | Line 96 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
96
97		}
98
99	<
100	<	//retrive output filename of z force
94	<	data = info->getProperty("zconsfilename");
99	>	//retrieve output filename of z force
100	>	data = info->getProperty(ZCONSFILENAME_ID);
101		if(!data) {
102
103
#	Line 104 | Line 110 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
110		}
111		else{
112
113	<	filename = dynamic_cast<StringData*>(data);
113	>	filename = dynamic_cast<StringData*>(data);
114
115		if(!filename){
116
#	Line 118 | Line 124 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
124		this->zconsOutput = filename->getData();
125		}
126
121	–
127		}
128
129	<
130	<	//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];
129	>	//retrieve tolerance for z-constraint molecuels
130	>	data = info->getProperty(ZCONSTOL_ID);
131
132	<	totalMass_local = 0;
133	<	totalMass = 0;
134	<	totalMZ_local = 0;
135	<	totalMZ = 0;
136	<	massOfUncons_local = 0;
137	<
132	>	if(!data) {
133	>
134	>	sprintf( painCave.errMsg,
135	>	"ZConstraint error: can not get tolerance \n");
136	>	painCave.isFatal = 1;
137	>	simError();
138	>	}
139	>	else{
140
141	<	for(int i = 0; i < nMols; i++){
142	<	massOfCurMol = molecules[i].getTotalMass();
143	<	molecules[i].getCOM(COM);
141	>	tolerance = dynamic_cast<DoubleData*>(data);
142
143	<	totalMass_local += massOfCurMol;
144	<	totalMZ_local += massOfCurMol * COM[2];
145	<
146	<	if(isZConstraintMol(&molecules[i]) == -1){
147	<
148	<	massOfUncons_local += massOfCurMol;
149	<	}
150	<
143	>	if(!tolerance){
144	>
145	>	sprintf( painCave.errMsg,
146	>	"ZConstraint error: Can not get property from SimInfo\n");
147	>	painCave.isFatal = 1;
148	>	simError();
149	>
150	>	}
151	>	else{
152	>	this->zconsTol = tolerance->getData();
153	>	}
154	>
155		}
156
157	<
158	<	#ifdef IS_MPI
159	<	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
157	>	//retrieve index of z-constraint molecules
158	>	data = info->getProperty(ZCONSPARADATA_ID);
159	>	if(!data) {
160
161	<	double zsys;
162	<	zsys = totalMZ / totalMass;
163	<
164	<	#ifndef IS_MPI
165	<	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	<
161	>	sprintf( painCave.errMsg,
162	>	"ZConstraint error: If you use an ZConstraint\n"
163	>	" , you must set index of z-constraint molecules.\n");
164	>	painCave.isFatal = 1;
165	>	simError();
166		}
167	<	#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;
167	>	else{
168
169	<	if(worldRank == 0){
169	>	zConsParaData = dynamic_cast<ZConsParaData*>(data);
170
171	<	int globalIndexOfCurMol;
172	<	int *MolToProcMap;
173	<	MolToProcMap = mpiSim->getMolToProcMap();
171	>	if(!zConsParaData){
172	>
173	>	sprintf( painCave.errMsg,
174	>	"ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
175	>	painCave.isFatal = 1;
176	>	simError();
177
178	<	for(int i = 0; i < indexOfAllZConsMols.size(); i++){
179	<
180	<	whichNode = MolToProcMap[indexOfAllZConsMols[i]];
181	<	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;
178	>	}
179	>	else{
180	>
181	>	parameters = zConsParaData->getData();
182
183	<	while (!done){
183	>	//check the range of zconsIndex
184	>	//and the minimum value of index is the first one (we already sorted the data)
185	>	//the maximum value of index is the last one
186	>
187	>	int maxIndex;
188	>	int minIndex;
189	>	int totalNumMol;
190	>
191	>	minIndex = (*parameters)[0].zconsIndex;
192	>	if(minIndex < 0){
193	>	sprintf( painCave.errMsg,
194	>	"ZConstraint error: index is out of range\n");
195	>	painCave.isFatal = 1;
196	>	simError();
197	>	}
198	>
199	>	maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
200	>
201	>	#ifndef IS_MPI
202	>	totalNumMol = nMols;
203	>	#else
204	>	totalNumMol = mpiSim->getTotNmol();
205	>	#endif
206
207	<	MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
208	<
209	<	switch (status){
210	<
211	<	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;
207	>	if(maxIndex > totalNumMol - 1){
208	>	sprintf( painCave.errMsg,
209	>	"ZConstraint error: index is out of range\n");
210	>	painCave.isFatal = 1;
211	>	simError();
212		}
261	–
262	–	}
263	–
264	–	}
213
214	<	//Brocast the allRefZ to slave nodes;
215	<	double* allRefZBuf;
216	<	int nZConsMols;
269	<	nZConsMols = indexOfAllZConsMols.size();
270	<
271	<	allRefZBuf = new double[nZConsMols];
272	<
273	<	if(worldRank == 0){
214	>	//if user does not specify the zpos for the zconstraint molecule
215	>	//its initial z coordinate  will be used as default
216	>	for(int i = 0; i < parameters->size(); i++){
217
218	<	for(int i = 0; i < nZConsMols; i++)
219	<	allRefZBuf[i] = allRefZ[i];
220	<	}
221	<
222	<	MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
223	<
224	<	if(worldRank != 0){
225	<
226	<	for(int i = 0; i < nZConsMols; i++)
227	<	allRefZ.push_back(allRefZBuf[i]);
228	<	}
229	<
230	<	delete[] allRefZBuf;
218	>	if(!(*parameters)[i].havingZPos){
219	>	#ifndef IS_MPI
220	>	for(int j = 0; j < nMols; j++){
221	>	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
222	>	molecules[j].getCOM(COM);
223	>	break;
224	>	}
225	>	}
226	>	#else
227	>	//query which processor current zconstraint molecule belongs to
228	>	int *MolToProcMap;
229	>	int whichNode;
230	>	double initZPos;
231	>	MolToProcMap = mpiSim->getMolToProcMap();
232	>	whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
233	>
234	>	//broadcast the zpos of current z-contraint molecule
235	>	//the node which contain this
236	>
237	>	if (worldRank == whichNode ){
238	>
239	>	for(int j = 0; j < nMols; j++)
240	>	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
241	>	molecules[j].getCOM(COM);
242	>	break;
243	>	}
244	>
245	>	}
246	>
247	>	MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);
248		#endif
249	+
250	+	(*parameters)[i].zPos = COM[whichDirection];
251
252	<
252	>	sprintf( painCave.errMsg,
253	>	"ZConstraint warning: Does not specify zpos for z-constraint molecule "
254	>	"initial z coornidate will be used \n");
255	>	painCave.isFatal = 0;
256	>	simError();
257	>
258	>	}
259	>	}
260	>
261	>	}//end if (!zConsParaData)
262	>	}//end  if (!data)
263	>
264	>	//
265		#ifdef IS_MPI
266		update();
267		#else
268		int searchResult;
269	<
296	<	refZ = allRefZ;
297	<
269	>
270		for(int i = 0; i < nMols; i++){
271
272		searchResult = isZConstraintMol(&molecules[i]);
#	Line 303 | Line 275 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
275
276		zconsMols.push_back(&molecules[i]);
277		massOfZConsMols.push_back(molecules[i].getTotalMass());
278	<
278	>
279	>	zPos.push_back((*parameters)[searchResult].zPos);
280	>	//       cout << "index: "<< (*parameters)[searchResult].zconsIndex
281	>	//              <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
282	>	kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
283	>
284		molecules[i].getCOM(COM);
285		}
286		else
#	Line 316 | Line 293 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
293		}
294
295		fz = new double[zconsMols.size()];
296	+	curZPos = new double[zconsMols.size()];
297		indexOfZConsMols = new int [zconsMols.size()];
298
299	<	if(!fz || !indexOfZConsMols){
299	>	if(!fz || !curZPos || !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++)
306	>	//determine the states of z-constraint molecules
307	>	for(int i = 0; i < zconsMols.size(); i++){
308		indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
309	+
310	+	zconsMols[i]->getCOM(COM);
311	+	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
312	+	states.push_back(zcsFixed);
313	+	else
314	+	states.push_back(zcsMoving);
315	+	}
316
317		#endif
318	+
319	+	//get total masss of unconstraint molecules
320	+	double totalMassOfUncons_local;
321	+	totalMassOfUncons_local = 0;
322
323	<	fzOut = new ZConsWriter(zconsOutput.c_str());
323	>	for(int i = 0; i < unconsMols.size(); i++)
324	>	totalMassOfUncons_local += unconsMols[i]->getTotalMass();
325	>
326	>	#ifndef IS_MPI
327	>	totalMassOfUncons = totalMassOfUncons_local;
328	>	#else
329	>	MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1,
330	>	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
331	>	#endif
332	>
333	>
334	>	//get total number of unconstrained atoms
335	>	int nUnconsAtoms_local;
336	>	nUnconsAtoms_local = 0;
337	>	for(int i = 0; i < unconsMols.size(); i++)
338	>	nUnconsAtoms_local += unconsMols[i]->getNAtoms();
339	>
340	>	#ifndef IS_MPI
341	>	totNumOfUnconsAtoms = nUnconsAtoms_local;
342	>	#else
343	>	MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1,
344	>	MPI_INT,MPI_SUM, MPI_COMM_WORLD);
345	>	#endif
346	>
347	>	// creat zconsWriter
348	>	fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);
349
350		if(!fzOut){
351		sprintf( painCave.errMsg,
#	Line 338 | Line 353 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
353		painCave.isFatal = 1;
354		simError();
355		}
356	<
357	<	fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
356	>
357	>	forcePolicy->update();
358		}
359
360		template<typename T> ZConstraint<T>::~ZConstraint()
361		{
362		if(fz)
363		delete[] fz;
364	+
365	+	if(curZPos)
366	+	delete[] curZPos;
367
368		if(indexOfZConsMols)
369		delete[] indexOfZConsMols;
370
371		if(fzOut)
372		delete fzOut;
373	+
374	+	if(forcePolicy)
375	+	delete forcePolicy;
376		}
377
378	+
379	+	/**
380	+	*
381	+	*/
382	+
383		#ifdef IS_MPI
384		template<typename T> void ZConstraint<T>::update()
385		{
#	Line 362 | Line 388 | template<typename T> void ZConstraint<T>::update()
388
389		zconsMols.clear();
390		massOfZConsMols.clear();
391	<	refZ.clear();
391	>	zPos.clear();
392	>	kz.clear();
393
394		unconsMols.clear();
395		massOfUnconsMols.clear();
#	Line 376 | Line 403 | template<typename T> void ZConstraint<T>::update()
403		if(index > -1){
404
405		zconsMols.push_back(&molecules[i]);
406	+	zPos.push_back((*parameters)[index].zPos);
407	+	kz.push_back((*parameters)[index].kRatio * zForceConst);
408	+
409		massOfZConsMols.push_back(molecules[i].getTotalMass());
410
411		molecules[i].getCOM(COM);
382	–	refZ.push_back(allRefZ[index]);
412		}
413		else
414		{
#	Line 389 | Line 418 | template<typename T> void ZConstraint<T>::update()
418
419		}
420		}
421	+
422	+	//determine the states of z-constraint molecules
423	+	for(int i = 0; i < zconsMols.size(); i++){
424	+	zconsMols[i]->getCOM(COM);
425	+	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
426	+	states.push_back(zcsFixed);
427	+	else
428	+	states.push_back(zcsMoving);
429	+	}
430	+
431
432		//The reason to declare fz and indexOfZconsMols as pointer to array is
433		// that we want to make the MPI communication simple
434		if(fz)
435		delete[] fz;
436	+
437	+	if(curZPos)
438	+	delete[] curZPos;
439
440		if(indexOfZConsMols)
441		delete[] indexOfZConsMols;
442
443		if (zconsMols.size() > 0){
444		fz = new double[zconsMols.size()];
445	+	curZPos = new double[zconsMols.size()];
446		indexOfZConsMols =  new int[zconsMols.size()];
447
448	<	if(!fz || !indexOfZConsMols){
448	>	if(!fz || !curZPos || !indexOfZConsMols){
449		sprintf( painCave.errMsg,
450		"Memory allocation failure in class Zconstraint\n");
451		painCave.isFatal = 1;
#	Line 416 | Line 459 | template<typename T> void ZConstraint<T>::update()
459		}
460		else{
461		fz = NULL;
462	+	curZPos = NULL;
463		indexOfZConsMols = NULL;
464		}
465
466	+	//
467	+	forcePolicy->update();
468	+
469		}
470
471		#endif
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
473	>	/**
474	>	*  Function Name: isZConstraintMol
475	>	*  Parameter
476	>	*    Molecule* mol
477	>	*  Return value:
478	>	*    -1, if the molecule is not z-constraint molecule,
479	>	*    other non-negative values, its index in indexOfAllZConsMols vector
480		*/
481
482		template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
#	Line 441 | Line 489 | template<typename T> int ZConstraint<T>::isZConstraint
489		index = mol->getGlobalIndex();
490
491		low = 0;
492	<	high = indexOfAllZConsMols.size() - 1;
492	>	high = parameters->size() - 1;
493
494		//Binary Search (we have sorted the array)
495		while(low <= high){
496		mid = (low + high) /2;
497	<	if (indexOfAllZConsMols[mid] == index)
497	>	if ((*parameters)[mid].zconsIndex == index)
498		return mid;
499	<	else if (indexOfAllZConsMols[mid] > index )
499	>	else if ((*parameters)[mid].zconsIndex > index )
500		high = mid -1;
501		else
502		low = mid + 1;
#	Line 457 | Line 505 | template<typename T> int ZConstraint<T>::isZConstraint
505		return -1;
506		}
507
508	<	/** 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();
508	>	template<typename T> void ZConstraint<T>::integrate(){
509
510	<	double currZConsTime = 0;
511	<
512	<	//write out forces of z constraint
482	<	if( info->getTime() >= currZConsTime){
483	<	fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
484	<	}
485	<	}
510	>	//zero out the velocities of center of mass of unconstrained molecules
511	>	//and the velocities of center of mass of every single z-constrained molecueles
512	>	zeroOutVel();
513
514	<	/** Function Name: resetZ
515	<	** Description:
516	<	**  Reset the z coordinates
490	<	*/
514	>	curZconsTime = zconsTime + info->getTime();
515	>
516	>	T::integrate();
517
518	<	template<typename T> void ZConstraint<T>::resetZ()
519	<	{
518	>	}
519	>
520
521	<	double pos[3];
522	<	double deltaZ;
523	<	double mzOfZCons;   //total sum of m*z of z-constrain molecules
524	<	double mzOfUncons; //total sum of m*z of unconstrain molecuels;
525	<	double totalMZOfZCons;
526	<	double totalMZOfUncons;
527	<	double COM[3];
521	>	/**
522	>	*
523	>	*
524	>	*
525	>	*
526	>	*/
527	>	template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
528		double zsys;
529	<	Atom** zconsAtoms;
529	>	double COM[3];
530	>	double force[3];
531	>	double zSysCOMVel;
532
533	<	mzOfZCons = 0;
534	<	mzOfUncons  = 0;
533	>	T::calcForce(calcPot, calcStress);
534	>
535	>	if (checkZConsState()){
536	>
537	>	#ifdef IS_MPI
538	>	if(worldRank == 0){
539	>	#endif
540	>	//       std::cerr << "\n"
541	>	//              << "*******************************************\n"
542	>	//              << " about to call zeroOutVel()\n"
543	>	//              << "*******************************************\n"
544	>	//              << "\n";
545	>	#ifdef IS_MPI
546	>	}
547	>	#endif
548	>	zeroOutVel();
549	>
550	>	#ifdef IS_MPI
551	>	if(worldRank == 0){
552	>	#endif
553	>	//       std::cerr << "\n"
554	>	//              << "*******************************************\n"
555	>	//              << " finished zeroOutVel()\n"
556	>	//              << "*******************************************\n"
557	>	//              << "\n";
558	>	#ifdef IS_MPI
559	>	}
560	>	#endif
561	>
562	>	forcePolicy->update();
563	>	}
564
565	<	for(int i = 0; i < zconsMols.size(); i++){
566	<	mzOfZCons += massOfZConsMols[i] * refZ[i];
567	<	}
565	>	zsys = calcZSys();
566	>	zSysCOMVel = calcSysCOMVel();
567	>	#ifdef IS_MPI
568	>	if(worldRank == 0){
569	>	#endif
570	>	//     cout << "---------------------------------------------------------------------" <<endl;
571	>	//     cout << "current time: " << info->getTime() << endl;
572	>	//     cout << "center of mass at z: " << zsys << endl;
573	>	//     cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
574
575		#ifdef IS_MPI
576	<	MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
514	<	#else
515	<	totalMZOfZCons = mzOfZCons;
576	>	}
577		#endif
578
579	<	for(int i = 0; i < unconsMols.size(); i++){
580	<	unconsMols[i]->getCOM(COM);
581	<	mzOfUncons += massOfUnconsMols[i] * COM[2];
579	>	//do zconstraint force;
580	>	if (haveFixedZMols())
581	>	this->doZconstraintForce();
582	>
583	>	//use harmonical poteintial to move the molecules to the specified positions
584	>	if (haveMovingZMols())
585	>	this->doHarmonic();
586	>
587	>	//write out forces and current positions of z-constraint molecules
588	>	if(info->getTime() >= curZconsTime){
589	>	for(int i = 0; i < zconsMols.size(); i++){
590	>	zconsMols[i]->getCOM(COM);
591	>	curZPos[i] = COM[whichDirection];
592	>
593	>	//if the z-constraint molecule is still moving, just record its force
594	>	if(states[i] == zcsMoving){
595	>	fz[i] = 0;
596	>	Atom** movingZAtoms;
597	>	movingZAtoms = zconsMols[i]->getMyAtoms();
598	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
599	>	movingZAtoms[j]->getFrc(force);
600	>	fz[i] += force[whichDirection];
601	>	}
602	>	}
603	>	}
604	>	fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos);
605	>	curZconsTime += zconsTime;
606		}
607	<
607	>
608	>	zSysCOMVel = calcSysCOMVel();
609		#ifdef IS_MPI
610	<	MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
611	<	#else
612	<	totalMZOfUncons = mzOfUncons;
613	<	#endif
614	<
615	<	zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons;
610	>	if(worldRank == 0){
611	>	#endif
612	>	//    cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
613	>	#ifdef IS_MPI
614	>	}
615	>	#endif
616	>	}
617	>
618	>
619	>	/**
620	>	*
621	>	*/
622
623	<	for(int i = 0; i < zconsMols.size(); i++){
623	>	template<typename T> double ZConstraint<T>::calcZSys()
624	>	{
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	>
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	>	}
644	>
645	>
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	>	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	>	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++){
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	<	zconsMols[i]->getCOM(COM);
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	<	deltaZ = zsys + refZ[i] - COM[2];
777	<	//update z coordinate
778	<	zconsAtoms = zconsMols[i]->getMyAtoms();
779	<	for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
780	<	zconsAtoms[j]->getPos(pos);
781	<	pos[2] += deltaZ;
782	<	zconsAtoms[j]->setPos(pos);
783	<	}
776	>	}
777	>
778	>	}
779	>
780	>
781	>	zSysCOMVel = calcSysCOMVel();
782	>	#ifdef IS_MPI
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	>	//constrain the molecules which do not reach the specified positions
806
807	<	//calculate z constrain force
808	<	fz[i] = massOfZConsMols[i]* deltaZ / dt2;
807	>	//Zero Out the force of z-contrained molecules
808	>	totalFZ_local = 0;
809	>
810	>	//calculate the total z-contrained force of fixed z-contrained molecules
811	>
812	>	for(int i = 0; i < zconsMols.size(); i++){
813	>
814	>	if (states[i] == zcsFixed){
815	>
816	>	zconsMols[i]->getCOM(COM);
817	>	zconsAtoms = zconsMols[i]->getMyAtoms();
818	>
819	>	fz[i] = 0;
820	>	for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
821	>	zconsAtoms[j]->getFrc(force);
822	>	fz[i] += force[whichDirection];
823	>	}
824	>	totalFZ_local += fz[i];
825	>
826	>	//cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
827	>	//      <<"\tcurrent zpos: " << COM[whichDirection]
828	>	//       << "\tcurrent fz: " <<fz[i] << endl;
829	>
830	>
831	>	}
832	>
833	>	}
834	>
835	>	//calculate total z-constraint force
836	>	#ifdef IS_MPI
837	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
838	>	#else
839	>	totalFZ = totalFZ_local;
840	>	#endif
841	>
842	>
843	>	// apply negative to fixed z-constrained molecues;
844	>	force[0]= 0;
845	>	force[1]= 0;
846	>	force[2]= 0;
847	>
848	>	for(int i = 0; i < zconsMols.size(); i++){
849	>
850	>	if (states[i] == zcsFixed){
851	>
852	>	int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
853	>	zconsAtoms = zconsMols[i]->getMyAtoms();
854	>
855	>	for(int j =0; j < nAtomOfCurZConsMol; j++) {
856	>	//force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
857	>	force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
858	>	zconsAtoms[j]->addFrc(force);
859	>	}
860	>
861	>	}
862	>
863	>	}
864	>
865	>	//   cout << "after zero out z-constraint force on fixed z-constraint molecuels "
866	>	//        << "total force is " << calcTotalForce() << endl;
867	>
868	>	force[0]= 0;
869	>	force[1]= 0;
870	>	force[2]= 0;
871	>
872	>	//modify the forces of unconstrained molecules
873	>	for(int i = 0; i < unconsMols.size(); i++){
874	>
875	>	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
876	>
877	>	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
878	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
879	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
880	>	unconsAtoms[j]->addFrc(force);
881	>	}
882	>
883	>	}
884	>
885	>	//modify the forces of moving z-constrained molecules
886	>	for(int i = 0; i < zconsMols.size(); i++) {
887	>	if (states[i] == zcsMoving){
888	>
889	>	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
890	>
891	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
892	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
893	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
894	>	movingZAtoms[j]->addFrc(force);
895	>	}
896	>	}
897	>	}
898	>
899	>	//cout << "after substracting z-constraint force from moving molecuels "
900	>	//      << "total force is " << calcTotalForce()  << endl;
901	>
902	>	}
903	>
904	>	/**
905	>	*
906	>	*
907	>	*/
908	>
909	>	template<typename T> void ZConstraint<T>::doHarmonic(){
910	>	double force[3];
911	>	double harmonicU;
912	>	double harmonicF;
913	>	double COM[3];
914	>	double diff;
915	>	double totalFZ_local;
916	>	double totalFZ;
917	>
918	>	force[0] = 0;
919	>	force[1] = 0;
920	>	force[2] = 0;
921	>
922	>	totalFZ_local = 0;
923	>
924	>	for(int i = 0; i < zconsMols.size(); i++) {
925	>
926	>	if (states[i] == zcsMoving){
927	>	zconsMols[i]->getCOM(COM);
928	>	//       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
929	>	//         << "\tcurrent zpos: " << COM[whichDirection] << endl;
930	>
931	>	diff = COM[whichDirection] -zPos[i];
932	>
933	>	harmonicU = 0.5 * kz[i] * diff * diff;
934	>	info->lrPot += harmonicU;
935	>
936	>	harmonicF =  - kz[i] * diff;
937	>	totalFZ_local += harmonicF;
938	>
939	>	//adjust force
940
941	+	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
942	+
943	+	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
944	+	//force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
945	+	force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
946	+	movingZAtoms[j]->addFrc(force);
947	+	}
948	+	}
949	+
950		}
951
952	+	#ifndef IS_MPI
953	+	totalFZ = totalFZ_local;
954	+	#else
955	+	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
956	+	#endif
957	+
958	+	force[0]= 0;
959	+	force[1]= 0;
960	+	force[2]= 0;
961	+
962	+	//modify the forces of unconstrained molecules
963	+	for(int i = 0; i < unconsMols.size(); i++){
964	+
965	+	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
966	+
967	+	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
968	+	//force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
969	+	force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
970	+	unconsAtoms[j]->addFrc(force);
971	+	}
972	+	}
973	+
974	+	}
975	+
976	+	/**
977	+	*
978	+	*/
979	+
980	+	template<typename T> bool ZConstraint<T>::checkZConsState(){
981	+	double COM[3];
982	+	double diff;
983	+
984	+	int changed_local;
985	+	int changed;
986	+
987	+	changed_local = 0;
988	+
989	+	for(int i =0; i < zconsMols.size(); i++){
990	+
991	+	zconsMols[i]->getCOM(COM);
992	+	diff = fabs(COM[whichDirection] - zPos[i]);
993	+	if (  diff <= zconsTol && states[i] == zcsMoving){
994	+	states[i] = zcsFixed;
995	+	changed_local = 1;
996	+	}
997	+	else if ( diff > zconsTol && states[i] == zcsFixed){
998	+	states[i] = zcsMoving;
999	+	changed_local = 1;
1000	+	}
1001	+
1002	+	}
1003	+
1004	+	#ifndef IS_MPI
1005	+	changed =changed_local;
1006	+	#else
1007	+	MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1008	+	#endif
1009	+
1010	+	return (changed > 0);
1011	+
1012	+	}
1013	+
1014	+	template<typename T> bool ZConstraint<T>::haveFixedZMols(){
1015	+
1016	+	int havingFixed_local;
1017	+	int havingFixed;
1018	+
1019	+	havingFixed_local = 0;
1020	+
1021	+	for(int i = 0; i < zconsMols.size(); i++)
1022	+	if (states[i] == zcsFixed){
1023	+	havingFixed_local = 1;
1024	+	break;
1025	+	}
1026	+
1027	+	#ifndef IS_MPI
1028	+	havingFixed = havingFixed_local;
1029	+	#else
1030	+	MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1031	+	#endif
1032	+
1033	+	return (havingFixed > 0);
1034	+	}
1035	+
1036	+
1037	+	/**
1038	+	*
1039	+	*/
1040	+	template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1041	+
1042	+	int havingMoving_local;
1043	+	int havingMoving;
1044	+
1045	+	havingMoving_local = 0;
1046	+
1047	+	for(int i = 0; i < zconsMols.size(); i++)
1048	+	if (states[i] == zcsMoving){
1049	+	havingMoving_local = 1;
1050	+	break;
1051	+	}
1052	+
1053	+	#ifndef IS_MPI
1054	+	havingMoving = havingMoving_local;
1055	+	#else
1056	+	MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1057	+	#endif
1058	+
1059	+	return (havingMoving > 0);
1060	+
1061	+	}
1062	+
1063	+	/**
1064	+	*
1065	+	*/
1066	+
1067	+	template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
1068	+	{
1069	+	double MVzOfMovingMols_local;
1070	+	double MVzOfMovingMols;
1071	+	double totalMassOfMovingZMols_local;
1072	+	double totalMassOfMovingZMols;
1073	+	double COMvel[3];
1074
1075	+	MVzOfMovingMols_local = 0;
1076	+	totalMassOfMovingZMols_local = 0;
1077	+
1078	+	for(int i =0; i < unconsMols.size(); i++){
1079	+	unconsMols[i]->getCOMvel(COMvel);
1080	+	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1081	+	}
1082	+
1083	+	for(int i = 0; i < zconsMols.size(); i++){
1084	+
1085	+	if (states[i] == zcsMoving){
1086	+	zconsMols[i]->getCOMvel(COMvel);
1087	+	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
1088	+	totalMassOfMovingZMols_local += massOfZConsMols[i];
1089	+	}
1090	+
1091	+	}
1092	+
1093	+	#ifndef IS_MPI
1094	+	MVzOfMovingMols = MVzOfMovingMols_local;
1095	+	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1096	+	#else
1097	+	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1098	+	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1099	+	#endif
1100	+
1101	+	double vzOfMovingMols;
1102	+	vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
1103	+
1104	+	return vzOfMovingMols;
1105		}
1106	+
1107	+	/**
1108	+	*
1109	+	*/
1110	+
1111	+	template<typename T> double ZConstraint<T>::calcSysCOMVel()
1112	+	{
1113	+	double COMvel[3];
1114	+	double tempMVz_local;
1115	+	double tempMVz;
1116	+	double massOfZCons_local;
1117	+	double massOfZCons;
1118	+
1119	+
1120	+	tempMVz_local = 0;
1121	+
1122	+	for(int i =0 ; i < nMols; i++){
1123	+	molecules[i].getCOMvel(COMvel);
1124	+	tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
1125	+	}
1126	+
1127	+	massOfZCons_local = 0;
1128	+
1129	+	for(int i = 0; i < massOfZConsMols.size(); i++){
1130	+	massOfZCons_local += massOfZConsMols[i];
1131	+	}
1132	+	#ifndef IS_MPI
1133	+	massOfZCons = massOfZCons_local;
1134	+	tempMVz = tempMVz_local;
1135	+	#else
1136	+	MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1137	+	MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1138	+	#endif
1139	+
1140	+	return tempMVz /(totalMassOfUncons + massOfZCons);
1141	+	}
1142	+
1143	+	/**
1144	+	*
1145	+	*/
1146	+
1147	+	template<typename T> double ZConstraint<T>::calcTotalForce(){
1148	+
1149	+	double force[3];
1150	+	double totalForce_local;
1151	+	double totalForce;
1152	+
1153	+	totalForce_local = 0;
1154	+
1155	+	for(int i = 0; i < nAtoms; i++){
1156	+	atoms[i]->getFrc(force);
1157	+	totalForce_local += force[whichDirection];
1158	+	}
1159	+
1160	+	#ifndef IS_MPI
1161	+	totalForce = totalForce_local;
1162	+	#else
1163	+	MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1164	+	#endif
1165	+
1166	+	return totalForce;
1167	+
1168	+	}
1169	+
1170	+	/**
1171	+	*
1172	+	*/
1173	+
1174	+	template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1175	+	//calculate the number of atoms of moving z-constrained molecules
1176	+	int nMovingZAtoms_local;
1177	+	int nMovingZAtoms;
1178	+
1179	+	nMovingZAtoms_local = 0;
1180	+	for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1181	+	if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1182	+	nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1183	+
1184	+	#ifdef IS_MPI
1185	+	MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1186	+	#else
1187	+	nMovingZAtoms = nMovingZAtoms_local;
1188	+	#endif
1189	+	totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1190	+
1191	+	#ifdef IS_MPI
1192	+	if(worldRank == 0){
1193	+	#endif
1194	+	//    std::cerr << "\n"
1195	+	//            << "*******************************************\n"
1196	+	//            << " fiished Policy by numbr()\n"
1197	+	//            << "*******************************************\n"
1198	+	//            << "\n";
1199	+	#ifdef IS_MPI
1200	+	}
1201	+	#endif
1202	+	}
1203	+
1204	+	template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1205	+	return totalForce / mol->getNAtoms();
1206	+	}
1207	+
1208	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, double totalForce){
1209	+	return totalForce / totNumOfMovingAtoms;
1210	+	}
1211	+
1212	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1213	+	return totalForce / mol->getNAtoms();
1214	+	}
1215	+
1216	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, double totalForce){
1217	+	return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1218	+	}
1219	+
1220	+	/**
1221	+	*
1222	+	*/
1223	+
1224	+	template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1225	+	//calculate the number of atoms of moving z-constrained molecules
1226	+	double massOfMovingZAtoms_local;
1227	+	double massOfMovingZAtoms;
1228	+
1229	+	massOfMovingZAtoms_local = 0;
1230	+	for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1231	+	if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1232	+	massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1233	+
1234	+	#ifdef IS_MPI
1235	+	MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1236	+	#else
1237	+	massOfMovingZAtoms = massOfMovingZAtoms_local;
1238	+	#endif
1239	+	totMassOfMovingAtoms = massOfMovingZAtoms_local + zconsIntegrator->totalMassOfUncons;
1240	+	}
1241	+
1242	+	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1243	+	return totalForce * atom->getMass() / mol->getTotalMass();
1244	+	}
1245	+
1246	+	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols( Atom* atom, double totalForce){
1247	+	return totalForce * atom->getMass() / totMassOfMovingAtoms;
1248	+	}
1249	+
1250	+	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1251	+	return totalForce * atom->getMass() / mol->getTotalMass();
1252	+	}
1253	+
1254	+	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, double totalForce){
1255	+	return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1256	+	}
1257	+

Diff Legend

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