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 658 by tim, Thu Jul 31 15:35:07 2003 UTC vs.
Revision 738 by tim, Tue Sep 2 14:30:12 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), fzOut(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 Subtraction 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 Subtraction policy, "
35	>	"PolicyByMass is used\n");
36	>	painCave.isFatal = 0;
37	>	simError();
38	>
39	>	forcePolicy = (ForceSubtractionPolicy*) 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 = (ForceSubtractionPolicy*) new PolicyByMass(this);
52		}
53		else{
54	<	indexOfAllZConsMols = index->getIndexData();
54	>	if(policy->getData() == "BYNUMBER")
55	>	forcePolicy = (ForceSubtractionPolicy*) new PolicyByNumber(this);
56	>	else if(policy->getData() == "BYMASS")
57	>	forcePolicy = (ForceSubtractionPolicy*) new PolicyByMass(this);
58	>	else{
59	>	sprintf( painCave.errMsg,
60	>	"ZConstraint Warning: unknown force Subtraction policy, "
61	>	"PolicyByMass is used\n");
62	>	painCave.isFatal = 0;
63	>	simError();
64	>	forcePolicy = (ForceSubtractionPolicy*) new PolicyByMass(this);
65	>	}
66		}
38	–
67		}
68
41	–	//retrive sample time of z-contraint
42	–	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 67 | Line 96 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
96
97		}
98
99	<
100	<	//retrive output filename of z force
72	<	data = info->getProperty("zconsfilename");
99	>	//retrieve output filename of z force
100	>	data = info->getProperty(ZCONSFILENAME_ID);
101		if(!data) {
102
103
#	Line 82 | 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 96 | Line 124 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
124		this->zconsOutput = filename->getData();
125		}
126
99	–
127		}
128
129	<
130	<	//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];
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	<
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	<
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	<
142	<	#ifdef IS_MPI
143	<	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
141	>	tolerance = dynamic_cast<DoubleData*>(data);
142	>
143	>	if(!tolerance){
144
145	<	double zsys;
146	<	zsys = totalMZ / totalMass;
147	<
148	<	#ifndef IS_MPI
149	<	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);
145	>	sprintf( painCave.errMsg,
146	>	"ZConstraint error: Can not get property from SimInfo\n");
147	>	painCave.isFatal = 1;
148	>	simError();
149	>
150		}
151	<
151	>	else{
152	>	this->zconsTol = tolerance->getData();
153	>	}
154	>
155		}
156	–	#else
157	–
158	–	int whichNode;
159	–	enum CommType { RequestMolZPos, EndOfRequest} status;
160	–	//int status;
161	–	double zpos;
162	–	int localIndex;
163	–	MPI_Status ierr;
164	–	int tag = 0;
156
157	<	if(worldRank == 0){
158	<
159	<	int globalIndexOfCurMol;
160	<	int *MolToProcMap;
161	<	MolToProcMap = mpiSim->getMolToProcMap();
162	<
163	<	for(int i = 0; i < indexOfAllZConsMols.size(); i++){
164	<
165	<	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	<
157	>	//retrieve index of z-constraint molecules
158	>	data = info->getProperty(ZCONSPARADATA_ID);
159	>	if(!data) {
160	>
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{
168
169	<	int whichMol;
170	<	bool done = false;
169	>	zConsParaData = dynamic_cast<ZConsParaData*>(data);
170	>
171	>	if(!zConsParaData){
172
173	<	while (!done){
174	<
175	<	MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
173	>	sprintf( painCave.errMsg,
174	>	"ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n");
175	>	painCave.isFatal = 1;
176	>	simError();
177
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	–
178		}
179	<
180	<	}
179	>	else{
180	>
181	>	parameters = zConsParaData->getData();
182
183	<	//Brocast the allRefZ to slave nodes;
184	<	double* allRefZBuf;
185	<	int nZConsMols;
247	<	nZConsMols = indexOfAllZConsMols.size();
248	<
249	<	allRefZBuf = new double[nZConsMols];
250	<
251	<	if(worldRank == 0){
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	<	for(int i = 0; i < nZConsMols; i++)
188	<	allRefZBuf[i] = allRefZ[i];
189	<	}
190	<
191	<	MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
192	<
193	<	if(worldRank != 0){
194	<
195	<	for(int i = 0; i < nZConsMols; i++)
196	<	allRefZ.push_back(allRefZBuf[i]);
197	<	}
198	<
199	<	delete[] allRefZBuf;
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	>	if(maxIndex > totalNumMol - 1){
208	>	sprintf( painCave.errMsg,
209	>	"ZConstraint error: index is out of range\n");
210	>	painCave.isFatal = 1;
211	>	simError();
212	>	}
213	>
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	>	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	<
274	<	refZ = allRefZ;
275	<
269	>
270		for(int i = 0; i < nMols; i++){
271
272		searchResult = isZConstraintMol(&molecules[i]);
#	Line 281 | 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 294 | 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());
324	<
325	<	if(!fzOut){
326	<	sprintf( painCave.errMsg,
327	<	"Memory allocation failure in class Zconstraint\n");
328	<	painCave.isFatal = 1;
329	<	simError();
330	<	}
331	<
332	<	fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
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	>	//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	>	forcePolicy->update();
347		}
348
349		template<typename T> ZConstraint<T>::~ZConstraint()
350		{
351		if(fz)
352		delete[] fz;
353	+
354	+	if(curZPos)
355	+	delete[] curZPos;
356
357		if(indexOfZConsMols)
358		delete[] indexOfZConsMols;
359
360		if(fzOut)
361		delete fzOut;
362	+
363	+	if(forcePolicy)
364	+	delete forcePolicy;
365		}
366
367	+
368	+	/**
369	+	*
370	+	*/
371	+
372		#ifdef IS_MPI
373		template<typename T> void ZConstraint<T>::update()
374		{
#	Line 340 | Line 377 | template<typename T> void ZConstraint<T>::update()
377
378		zconsMols.clear();
379		massOfZConsMols.clear();
380	<	refZ.clear();
380	>	zPos.clear();
381	>	kz.clear();
382
383		unconsMols.clear();
384		massOfUnconsMols.clear();
#	Line 354 | Line 392 | template<typename T> void ZConstraint<T>::update()
392		if(index > -1){
393
394		zconsMols.push_back(&molecules[i]);
395	+	zPos.push_back((*parameters)[index].zPos);
396	+	kz.push_back((*parameters)[index].kRatio * zForceConst);
397	+
398		massOfZConsMols.push_back(molecules[i].getTotalMass());
399
400		molecules[i].getCOM(COM);
360	–	refZ.push_back(allRefZ[index]);
401		}
402		else
403		{
#	Line 367 | Line 407 | template<typename T> void ZConstraint<T>::update()
407
408		}
409		}
410	+
411	+	//determine the states of z-constraint molecules
412	+	for(int i = 0; i < zconsMols.size(); i++){
413	+	zconsMols[i]->getCOM(COM);
414	+	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
415	+	states.push_back(zcsFixed);
416	+	else
417	+	states.push_back(zcsMoving);
418	+	}
419	+
420
421		//The reason to declare fz and indexOfZconsMols as pointer to array is
422		// that we want to make the MPI communication simple
423		if(fz)
424		delete[] fz;
425	+
426	+	if(curZPos)
427	+	delete[] curZPos;
428
429		if(indexOfZConsMols)
430		delete[] indexOfZConsMols;
431
432		if (zconsMols.size() > 0){
433		fz = new double[zconsMols.size()];
434	+	curZPos = new double[zconsMols.size()];
435		indexOfZConsMols =  new int[zconsMols.size()];
436
437	<	if(!fz || !indexOfZConsMols){
437	>	if(!fz || !curZPos || !indexOfZConsMols){
438		sprintf( painCave.errMsg,
439		"Memory allocation failure in class Zconstraint\n");
440		painCave.isFatal = 1;
#	Line 394 | Line 448 | template<typename T> void ZConstraint<T>::update()
448		}
449		else{
450		fz = NULL;
451	+	curZPos = NULL;
452		indexOfZConsMols = NULL;
453		}
454
455	+	//
456	+	forcePolicy->update();
457	+
458		}
459
460		#endif
461
462	<	/**  Function Name: isZConstraintMol
463	<	**  Parameter
464	<	**    Molecule* mol
465	<	**  Return value:
466	<	**    -1, if the molecule is not z-constraint molecule,
467	<	**    other non-negative values, its index in indexOfAllZConsMols vector
462	>	/**
463	>	*  Function Name: isZConstraintMol
464	>	*  Parameter
465	>	*    Molecule* mol
466	>	*  Return value:
467	>	*    -1, if the molecule is not z-constraint molecule,
468	>	*    other non-negative values, its index in indexOfAllZConsMols vector
469		*/
470
471		template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
#	Line 419 | Line 478 | template<typename T> int ZConstraint<T>::isZConstraint
478		index = mol->getGlobalIndex();
479
480		low = 0;
481	<	high = indexOfAllZConsMols.size() - 1;
481	>	high = parameters->size() - 1;
482
483		//Binary Search (we have sorted the array)
484		while(low <= high){
485		mid = (low + high) /2;
486	<	if (indexOfAllZConsMols[mid] == index)
486	>	if ((*parameters)[mid].zconsIndex == index)
487		return mid;
488	<	else if (indexOfAllZConsMols[mid] > index )
488	>	else if ((*parameters)[mid].zconsIndex > index )
489		high = mid -1;
490		else
491		low = mid + 1;
#	Line 435 | Line 494 | template<typename T> int ZConstraint<T>::isZConstraint
494		return -1;
495		}
496
497	<	/** Function Name: integrateStep
498	<	** Parameter:
499	<	**   int calcPot;
500	<	**   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();
497	>	template<typename T> void ZConstraint<T>::integrate(){
498	>
499	>	// creat zconsWriter
500	>	fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);
501
502	<	double currZConsTime = 0;
503	<
504	<	//write out forces of z constraint
505	<	if( info->getTime() >= currZConsTime){
506	<	fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
507	<	}
502	>	if(!fzOut){
503	>	sprintf( painCave.errMsg,
504	>	"Memory allocation failure in class Zconstraint\n");
505	>	painCave.isFatal = 1;
506	>	simError();
507	>	}
508	>
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	>	zeroOutVel();
536	>	forcePolicy->update();
537	>	}
538	>
539	>	zsys = calcZSys();
540	>	zSysCOMVel = calcSysCOMVel();
541	>	#ifdef IS_MPI
542	>	if(worldRank == 0){
543	>	#endif
544	>	//cout << "---------------------------------------------------------------------" <<endl;
545	>	//cout << "current time: " << info->getTime() << endl;
546	>	//cout << "center of mass at z: " << zsys << endl;
547	>	//cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
548	>
549	>	#ifdef IS_MPI
550	>	}
551	>	#endif
552	>
553	>	//do zconstraint force;
554	>	if (haveFixedZMols())
555	>	this->doZconstraintForce();
556	>
557	>	//use harmonical poteintial to move the molecules to the specified positions
558	>	if (haveMovingZMols())
559	>	this->doHarmonic();
560	>
561	>	//write out forces and current positions of z-constraint molecules
562	>	if(info->getTime() >= curZconsTime){
563	>	for(int i = 0; i < zconsMols.size(); i++){
564	>	zconsMols[i]->getCOM(COM);
565	>	curZPos[i] = COM[whichDirection];
566	>
567	>	//if the z-constraint molecule is still moving, just record its force
568	>	if(states[i] == zcsMoving){
569	>	fz[i] = 0;
570	>	Atom** movingZAtoms;
571	>	movingZAtoms = zconsMols[i]->getMyAtoms();
572	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
573	>	movingZAtoms[j]->getFrc(force);
574	>	fz[i] += force[whichDirection];
575	>	}
576	>	}
577	>	}
578	>	fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos);
579	>	curZconsTime += zconsTime;
580	>	}
581	>
582	>	zSysCOMVel = calcSysCOMVel();
583	>	#ifdef IS_MPI
584	>	if(worldRank == 0){
585	>	#endif
586	>	//cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
587	>	#ifdef IS_MPI
588	>	}
589	>	#endif
590	>	}
591	>
592	>
593	>	/**
594	>	*
595	>	*/
596	>
597	>	template<typename T> double ZConstraint<T>::calcZSys()
598		{
599	<	double deltaZ;
600	<	double mzOfZCons;   //total sum of m*z of z-constrain molecules
601	<	double mzOfUncons; //total sum of m*z of unconstrain molecuels;
602	<	double totalMZOfZCons;
603	<	double totalMZOfUncons;
599	>	//calculate reference z coordinate for z-constraint molecules
600	>	double totalMass_local;
601	>	double totalMass;
602	>	double totalMZ_local;
603	>	double totalMZ;
604	>	double massOfCurMol;
605		double COM[3];
606	+
607	+	totalMass_local = 0;
608	+	totalMZ_local = 0;
609	+
610	+	for(int i = 0; i < nMols; i++){
611	+	massOfCurMol = molecules[i].getTotalMass();
612	+	molecules[i].getCOM(COM);
613	+
614	+	totalMass_local += massOfCurMol;
615	+	totalMZ_local += massOfCurMol * COM[whichDirection];
616	+
617	+	}
618	+
619	+
620	+	#ifdef IS_MPI
621	+	MPI_Allreduce(&totalMass_local, &totalMass, 1,
622	+	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
623	+	MPI_Allreduce(&totalMZ_local, &totalMZ, 1,
624	+	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
625	+	#else
626	+	totalMass = totalMass_local;
627	+	totalMZ = totalMZ_local;
628	+	#endif
629	+
630		double zsys;
631	<	Atom** zconsAtoms;
631	>	zsys = totalMZ / totalMass;
632
633	<	mzOfZCons = 0;
634	<	mzOfUncons  = 0;
633	>	return zsys;
634	>	}
635	>
636	>	/**
637	>	*
638	>	*/
639	>	template<typename T> void ZConstraint<T>::thermalize( void ){
640	>
641	>	T::thermalize();
642	>	zeroOutVel();
643	>	}
644	>
645	>	/**
646	>	*
647	>	*/
648	>
649	>	template<typename T> void ZConstraint<T>::zeroOutVel(){
650	>
651	>	Atom** fixedZAtoms;
652	>	double COMvel[3];
653	>	double vel[3];
654	>	double zSysCOMVel;
655	>
656	>	//zero out the velocities of center of mass of fixed z-constrained molecules
657
658		for(int i = 0; i < zconsMols.size(); i++){
659	<	mzOfZCons += massOfZConsMols[i] * refZ[i];
659	>
660	>	if (states[i] == zcsFixed){
661	>
662	>	zconsMols[i]->getCOMvel(COMvel);
663	>	//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
664	>
665	>	fixedZAtoms = zconsMols[i]->getMyAtoms();
666	>
667	>	for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
668	>	fixedZAtoms[j]->getVel(vel);
669	>	vel[whichDirection] -= COMvel[whichDirection];
670	>	fixedZAtoms[j]->setVel(vel);
671	>	}
672	>
673	>	zconsMols[i]->getCOMvel(COMvel);
674	>	//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
675	>	}
676	>
677		}
678
679	+	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
680	+
681	+	zSysCOMVel = calcSysCOMVel();
682		#ifdef IS_MPI
683	<	MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
683	>	if(worldRank == 0){
684	>	#endif
685	>	//cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;
686	>	#ifdef IS_MPI
687	>	}
688	>	#endif
689	>
690	>	// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
691	>	double MVzOfMovingMols_local;
692	>	double MVzOfMovingMols;
693	>	double totalMassOfMovingZMols_local;
694	>	double totalMassOfMovingZMols;
695	>
696	>	MVzOfMovingMols_local = 0;
697	>	totalMassOfMovingZMols_local = 0;
698	>
699	>	for(int i =0; i < unconsMols.size(); i++){
700	>	unconsMols[i]->getCOMvel(COMvel);
701	>	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
702	>	}
703	>
704	>	for(int i = 0; i < zconsMols.size(); i++){
705	>	if (states[i] == zcsMoving){
706	>	zconsMols[i]->getCOMvel(COMvel);
707	>	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
708	>	totalMassOfMovingZMols_local += massOfZConsMols[i];
709	>	}
710	>
711	>	}
712	>
713	>	#ifndef IS_MPI
714	>	MVzOfMovingMols = MVzOfMovingMols_local;
715	>	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
716		#else
717	<	totalMZOfZCons = mzOfZCons;
717	>	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
718	>	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
719		#endif
720
721	+	double vzOfMovingMols;
722	+	vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
723	+
724	+	//modify the velocites of unconstrained molecules
725	+	Atom** unconsAtoms;
726		for(int i = 0; i < unconsMols.size(); i++){
727	<	unconsMols[i]->getCOM(COM);
728	<	mzOfUncons += massOfUnconsMols[i] * COM[2];
727	>
728	>	unconsAtoms = unconsMols[i]->getMyAtoms();
729	>	for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
730	>	unconsAtoms[j]->getVel(vel);
731	>	vel[whichDirection] -= vzOfMovingMols;
732	>	unconsAtoms[j]->setVel(vel);
733	>	}
734	>
735	>	}
736	>
737	>	//modify the velocities of moving z-constrained molecuels
738	>	Atom** movingZAtoms;
739	>	for(int i = 0; i < zconsMols.size(); i++){
740	>
741	>	if (states[i] ==zcsMoving){
742	>
743	>	movingZAtoms = zconsMols[i]->getMyAtoms();
744	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
745	>	movingZAtoms[j]->getVel(vel);
746	>	vel[whichDirection] -= vzOfMovingMols;
747	>	movingZAtoms[j]->setVel(vel);
748	>	}
749	>
750	>	}
751	>
752	>	}
753	>
754	>
755	>	zSysCOMVel = calcSysCOMVel();
756	>	#ifdef IS_MPI
757	>	if(worldRank == 0){
758	>	#endif
759	>	//cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;
760	>	#ifdef IS_MPI
761		}
762	+	#endif
763	+
764	+	}
765	+
766	+	/**
767	+	*
768	+	*/
769	+
770	+	template<typename T> void ZConstraint<T>::doZconstraintForce(){
771	+
772	+	Atom** zconsAtoms;
773	+	double totalFZ;
774	+	double totalFZ_local;
775	+	double COMvel[3];
776	+	double COM[3];
777	+	double force[3];
778	+
779	+	//constrain the molecules which do not reach the specified positions
780	+
781	+	//Zero Out the force of z-contrained molecules
782	+	totalFZ_local = 0;
783	+
784	+	//calculate the total z-contrained force of fixed z-contrained molecules
785
786	+	//cout << "before zero out z-constraint force on fixed z-constraint molecuels "
787	+	//       << "total force is " << calcTotalForce() << endl;
788	+
789	+	for(int i = 0; i < zconsMols.size(); i++){
790	+
791	+	if (states[i] == zcsFixed){
792	+
793	+	zconsMols[i]->getCOM(COM);
794	+	zconsAtoms = zconsMols[i]->getMyAtoms();
795	+
796	+	fz[i] = 0;
797	+	for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
798	+	zconsAtoms[j]->getFrc(force);
799	+	fz[i] += force[whichDirection];
800	+	}
801	+	totalFZ_local += fz[i];
802	+
803	+	//cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
804	+	//      <<"\tcurrent zpos: " << COM[whichDirection]
805	+	//      << "\tcurrent fz: " <<fz[i] << endl;
806	+
807	+
808	+	}
809	+
810	+	}
811	+
812	+	//calculate total z-constraint force
813		#ifdef IS_MPI
814	<	MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
814	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
815		#else
816	<	totalMZOfUncons = mzOfUncons;
817	<	#endif
816	>	totalFZ = totalFZ_local;
817	>	#endif
818	>
819
820	<	zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons;
820	>	// apply negative to fixed z-constrained molecues;
821	>	force[0]= 0;
822	>	force[1]= 0;
823	>	force[2]= 0;
824
825	<	cout << "current time: " << info->getTime() <<endl;
826	<	for(int i = 0; i < zconsMols.size(); i++){
825	>	for(int i = 0; i < zconsMols.size(); i++){
826	>
827	>	if (states[i] == zcsFixed){
828	>
829	>	int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
830	>	zconsAtoms = zconsMols[i]->getMyAtoms();
831	>
832	>	for(int j =0; j < nAtomOfCurZConsMol; j++) {
833	>	//force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
834	>	force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
835	>	zconsAtoms[j]->addFrc(force);
836	>	}
837	>
838	>	}
839	>
840	>	}
841	>
842	>	//cout << "after zero out z-constraint force on fixed z-constraint molecuels "
843	>	//      << "total force is " << calcTotalForce() << endl;
844
845	<	zconsMols[i]->getCOM(COM);
845	>
846	>	force[0]= 0;
847	>	force[1]= 0;
848	>	force[2]= 0;
849	>
850	>	//modify the forces of unconstrained molecules
851	>	for(int i = 0; i < unconsMols.size(); i++){
852	>
853	>	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
854	>
855	>	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
856	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
857	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
858	>	unconsAtoms[j]->addFrc(force);
859	>	}
860	>
861	>	}
862	>
863	>	//modify the forces of moving z-constrained molecules
864	>	for(int i = 0; i < zconsMols.size(); i++) {
865	>	if (states[i] == zcsMoving){
866
867	<	cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t";
868	<	deltaZ = zsys + refZ[i] - COM[2];
869	<	cout << "\tdistance: " << COM[2] +deltaZ - zsys;
870	<	//update z coordinate
871	<	zconsAtoms = zconsMols[i]->getMyAtoms();
872	<	for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
873	<	zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ);
874	<	}
867	>	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
868	>
869	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
870	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
871	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
872	>	movingZAtoms[j]->addFrc(force);
873	>	}
874	>	}
875	>	}
876	>	//  cout << "after substracting z-constraint force from moving molecuels "
877	>	//        << "total force is " << calcTotalForce()  << endl;
878	>
879	>	}
880	>
881	>	/**
882	>	*
883	>	*
884	>	*/
885	>
886	>	template<typename T> void ZConstraint<T>::doHarmonic(){
887	>	double force[3];
888	>	double harmonicU;
889	>	double harmonicF;
890	>	double COM[3];
891	>	double diff;
892	>	double totalFZ_local;
893	>	double totalFZ;
894	>
895	>	force[0] = 0;
896	>	force[1] = 0;
897	>	force[2] = 0;
898	>
899	>	totalFZ_local = 0;
900	>
901	>	for(int i = 0; i < zconsMols.size(); i++) {
902	>
903	>	if (states[i] == zcsMoving){
904	>	zconsMols[i]->getCOM(COM);
905	>	//       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
906	>	//         << "\tcurrent zpos: " << COM[whichDirection] << endl;
907	>
908	>	diff = COM[whichDirection] -zPos[i];
909
910	<	//calculate z constrain force
911	<	fz[i] = massOfZConsMols[i]* deltaZ / dt2;
910	>	harmonicU = 0.5 * kz[i] * diff * diff;
911	>	info->lrPot += harmonicU;
912	>
913	>	harmonicF =  - kz[i] * diff;
914	>	totalFZ_local += harmonicF;
915	>
916	>	//adjust force
917
918	<	cout << "\tforce: " << fz[i] << endl;
918	>	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
919	>
920	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
921	>	//force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
922	>	force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
923	>	movingZAtoms[j]->addFrc(force);
924	>	}
925	>	}
926	>
927		}
928
929	+	#ifndef IS_MPI
930	+	totalFZ = totalFZ_local;
931	+	#else
932	+	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
933	+	#endif
934	+
935	+	cout << "before substracting harmonic force from moving molecuels "
936	+	<< "total force is " << calcTotalForce()  << endl;
937	+
938	+	force[0]= 0;
939	+	force[1]= 0;
940	+	force[2]= 0;
941	+
942	+	//modify the forces of unconstrained molecules
943	+	for(int i = 0; i < unconsMols.size(); i++){
944	+
945	+	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
946	+
947	+	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
948	+	//force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
949	+	force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
950	+	unconsAtoms[j]->addFrc(force);
951	+	}
952	+	}
953	+
954	+	cout << "after substracting harmonic force from moving molecuels "
955	+	<< "total force is " << calcTotalForce()  << endl;
956	+
957	+	}
958	+
959	+	/**
960	+	*
961	+	*/
962	+
963	+	template<typename T> bool ZConstraint<T>::checkZConsState(){
964	+	double COM[3];
965	+	double diff;
966	+
967	+	int changed_local;
968	+	int changed;
969	+
970	+	changed_local = 0;
971	+
972	+	for(int i =0; i < zconsMols.size(); i++){
973	+
974	+	zconsMols[i]->getCOM(COM);
975	+	diff = fabs(COM[whichDirection] - zPos[i]);
976	+	if (  diff <= zconsTol && states[i] == zcsMoving){
977	+	states[i] = zcsFixed;
978	+	changed_local = 1;
979	+	}
980	+	else if ( diff > zconsTol && states[i] == zcsFixed){
981	+	states[i] = zcsMoving;
982	+	changed_local = 1;
983	+	}
984	+
985	+	}
986	+
987	+	#ifndef IS_MPI
988	+	changed =changed_local;
989	+	#else
990	+	MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
991	+	#endif
992	+
993	+	return (changed > 0);
994	+
995	+	}
996	+
997	+	template<typename T> bool ZConstraint<T>::haveFixedZMols(){
998	+
999	+	int havingFixed_local;
1000	+	int havingFixed;
1001	+
1002	+	havingFixed_local = 0;
1003	+
1004	+	for(int i = 0; i < zconsMols.size(); i++)
1005	+	if (states[i] == zcsFixed){
1006	+	havingFixed_local = 1;
1007	+	break;
1008	+	}
1009	+
1010	+	#ifndef IS_MPI
1011	+	havingFixed = havingFixed_local;
1012	+	#else
1013	+	MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1014	+	#endif
1015	+
1016	+	return (havingFixed > 0);
1017	+	}
1018	+
1019	+
1020	+	/**
1021	+	*
1022	+	*/
1023	+	template<typename T> bool ZConstraint<T>::haveMovingZMols(){
1024	+
1025	+	int havingMoving_local;
1026	+	int havingMoving;
1027	+
1028	+	havingMoving_local = 0;
1029	+
1030	+	for(int i = 0; i < zconsMols.size(); i++)
1031	+	if (states[i] == zcsMoving){
1032	+	havingMoving_local = 1;
1033	+	break;
1034	+	}
1035	+
1036	+	#ifndef IS_MPI
1037	+	havingMoving = havingMoving_local;
1038	+	#else
1039	+	MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1040	+	#endif
1041	+
1042	+	return (havingMoving > 0);
1043	+
1044	+	}
1045	+
1046	+	/**
1047	+	*
1048	+	*/
1049	+
1050	+	template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
1051	+	{
1052	+	double MVzOfMovingMols_local;
1053	+	double MVzOfMovingMols;
1054	+	double totalMassOfMovingZMols_local;
1055	+	double totalMassOfMovingZMols;
1056	+	double COMvel[3];
1057
1058	+	MVzOfMovingMols_local = 0;
1059	+	totalMassOfMovingZMols_local = 0;
1060	+
1061	+	for(int i =0; i < unconsMols.size(); i++){
1062	+	unconsMols[i]->getCOMvel(COMvel);
1063	+	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1064	+	}
1065	+
1066	+	for(int i = 0; i < zconsMols.size(); i++){
1067	+
1068	+	if (states[i] == zcsMoving){
1069	+	zconsMols[i]->getCOMvel(COMvel);
1070	+	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
1071	+	totalMassOfMovingZMols_local += massOfZConsMols[i];
1072	+	}
1073	+
1074	+	}
1075	+
1076	+	#ifndef IS_MPI
1077	+	MVzOfMovingMols = MVzOfMovingMols_local;
1078	+	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1079	+	#else
1080	+	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1081	+	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1082	+	#endif
1083	+
1084	+	double vzOfMovingMols;
1085	+	vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
1086	+
1087	+	return vzOfMovingMols;
1088		}
1089	+
1090	+	/**
1091	+	*
1092	+	*/
1093	+
1094	+	template<typename T> double ZConstraint<T>::calcSysCOMVel()
1095	+	{
1096	+	double COMvel[3];
1097	+	double tempMVz_local;
1098	+	double tempMVz;
1099	+	double massOfZCons_local;
1100	+	double massOfZCons;
1101	+
1102	+
1103	+	tempMVz_local = 0;
1104	+
1105	+	for(int i =0 ; i < nMols; i++){
1106	+	molecules[i].getCOMvel(COMvel);
1107	+	tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
1108	+	}
1109	+
1110	+	massOfZCons_local = 0;
1111	+
1112	+	for(int i = 0; i < massOfZConsMols.size(); i++){
1113	+	massOfZCons_local += massOfZConsMols[i];
1114	+	}
1115	+	#ifndef IS_MPI
1116	+	massOfZCons = massOfZCons_local;
1117	+	tempMVz = tempMVz_local;
1118	+	#else
1119	+	MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1120	+	MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1121	+	#endif
1122	+
1123	+	return tempMVz /(totalMassOfUncons + massOfZCons);
1124	+	}
1125	+
1126	+	/**
1127	+	*
1128	+	*/
1129	+
1130	+	template<typename T> double ZConstraint<T>::calcTotalForce(){
1131	+
1132	+	double force[3];
1133	+	double totalForce_local;
1134	+	double totalForce;
1135	+
1136	+	totalForce_local = 0;
1137	+
1138	+	for(int i = 0; i < nAtoms; i++){
1139	+	atoms[i]->getFrc(force);
1140	+	totalForce_local += force[whichDirection];
1141	+	}
1142	+
1143	+	#ifndef IS_MPI
1144	+	totalForce = totalForce_local;
1145	+	#else
1146	+	MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1147	+	#endif
1148	+
1149	+	return totalForce;
1150	+
1151	+	}
1152	+
1153	+	/**
1154	+	*
1155	+	*/
1156	+
1157	+	template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1158	+	//calculate the number of atoms of moving z-constrained molecules
1159	+	int nMovingZAtoms_local;
1160	+	int nMovingZAtoms;
1161	+
1162	+	nMovingZAtoms_local = 0;
1163	+	for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1164	+	if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1165	+	nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1166	+
1167	+	#ifdef IS_MPI
1168	+	MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1169	+	#else
1170	+	nMovingZAtoms = nMovingZAtoms_local;
1171	+	#endif
1172	+	totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1173	+	}
1174	+
1175	+	template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1176	+	return totalForce / mol->getNAtoms();
1177	+	}
1178	+
1179	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, double totalForce){
1180	+	return totalForce / totNumOfMovingAtoms;
1181	+	}
1182	+
1183	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1184	+	return totalForce / mol->getNAtoms();
1185	+	}
1186	+
1187	+	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, double totalForce){
1188	+	return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1189	+	}
1190	+
1191	+	/**
1192	+	*
1193	+	*/
1194	+
1195	+	template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1196	+	//calculate the number of atoms of moving z-constrained molecules
1197	+	double massOfMovingZAtoms_local;
1198	+	double massOfMovingZAtoms;
1199	+
1200	+	massOfMovingZAtoms_local = 0;
1201	+	for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1202	+	if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1203	+	massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1204	+
1205	+	#ifdef IS_MPI
1206	+	MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1207	+	#else
1208	+	massOfMovingZAtoms = massOfMovingZAtoms_local;
1209	+	#endif
1210	+	totMassOfMovingAtoms = massOfMovingZAtoms + zconsIntegrator->totalMassOfUncons;
1211	+	}
1212	+
1213	+	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1214	+	return totalForce * atom->getMass() / mol->getTotalMass();
1215	+	}
1216	+
1217	+	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols( Atom* atom, double totalForce){
1218	+	return totalForce * atom->getMass() / totMassOfMovingAtoms;
1219	+	}
1220	+
1221	+	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1222	+	return totalForce * atom->getMass() / mol->getTotalMass();
1223	+	}
1224	+
1225	+	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, double totalForce){
1226	+	return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1227	+	}
1228	+

Diff Legend

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