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

Diff Legend

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