[ViewVC] Diff of: group/trunk/OOPSE/libmdtools/ZConstraint.cpp

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 676 by tim, Mon Aug 11 19:40:06 2003 UTC vs.
Revision 830 by gezelter, Tue Oct 28 16:20:28 2003 UTC

#	Line 1 \| Line 1
1		#include "Integrator.hpp"
2		#include "simError.h"
3	<	#include <cmath>
3	>	#include <math.h>
4		template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff)
5	<	: T(theInfo, the_ff), fz(NULL),
6	<	indexOfZConsMols(NULL)
5	>	: T(theInfo, the_ff), indexOfZConsMols(NULL), fz(NULL), curZPos(NULL),
6	>	fzOut(NULL), curZconsTime(0), forcePolicy(NULL)
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	<	//retrieve index of z-constraint molecules
28	<	data = info->getProperty("zconsindex");
17	<	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){
29	<
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	<
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;
42	<	int minIndex;
43	<	int totalNumMol;
44	<
45	<	minIndex = indexOfAllZConsMols[0];
46	<	if(minIndex < 0){
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 error: index is out of range\n");
61	<	painCave.isFatal = 1;
62	<	simError();
63	<	}
64	<
65	<	maxIndex = indexOfAllZConsMols[indexOfAllZConsMols.size() - 1];
54	<
55	<	#ifndef IS_MPI
56	<	totalNumMol = nMols;
57	<	#else
58	<	totalNumMol = mpiSim->getTotNmol();
59	<	#endif
60	<
61	<	if(maxIndex > totalNumMol - 1){
62	<	sprintf( painCave.errMsg,
63	<	"ZConstraint error: index is out of range\n");
64	<	painCave.isFatal = 1;
65	<	simError();
66	<
67	<	}
68	<
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		}
70	–
67		}
68
69	+
70		//retrieve sample time of z-contraint
71	<	data = info->getProperty("zconstime");
71	>	data = info->getProperty(ZCONSTIME_ID);
72
73		if(!data) {
74
#	Line 99 \| Line 96 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
96
97		}
98
102	–
99		//retrieve output filename of z force
100	<	data = info->getProperty("zconsfilename");
100	>	data = info->getProperty(ZCONSFILENAME_ID);
101		if(!data) {
102
103
#	Line 114 \| 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 128 \| Line 124 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
124		this->zconsOutput = filename->getData();
125		}
126
127	+	}
128
129	+	//retrieve tolerance for z-constraint molecuels
130	+	data = info->getProperty(ZCONSTOL_ID);
131	+
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	+	tolerance = dynamic_cast<DoubleData*>(data);
142	+
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	+	//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	+	zConsParaData = dynamic_cast<ZConsParaData*>(data);
170	+
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	+	}
179	+	else{
180	+
181	+	parameters = zConsParaData->getData();
182	+
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	+	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 < (int)(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	+
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, whichNode, MPI_COMM_WORLD);
248	+	#endif
249	+
250	+	(*parameters)[i].zPos = COM[whichDirection];
251	+
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;
138	–	double COM[3];
269
270		for(int i = 0; i < nMols; i++){
271
#	Line 145 \| Line 275 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
275
276		zconsMols.push_back(&molecules[i]);
277		massOfZConsMols.push_back(molecules[i].getTotalMass());
278	+
279	+	zPos.push_back((*parameters)[searchResult].zPos);
280	+	// cout << "index: "<< (*parameters)[searchResult].zconsIndex
281	+	// <<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
282
283	+	kz.push_back((parameters)[searchResult]. kRatio zForceConst);
284		molecules[i].getCOM(COM);
285		}
286		else
#	Line 158 \| 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 < (int)(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	+	for(int i = 0; i < (int)(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++)
336	>	for(int i = 0; i < (int)(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, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
343	<	#endif
342	>	MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1,
343	>	MPI_INT,MPI_SUM, MPI_COMM_WORLD);
344	>	#endif
345
346	<
188	<
189	<	fzOut = new ZConsWriter(zconsOutput.c_str());
190	<
191	<	if(!fzOut){
192	<	sprintf( painCave.errMsg,
193	<	"Memory allocation failure in class Zconstraint\n");
194	<	painCave.isFatal = 1;
195	<	simError();
196	<	}
197	<
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 217 \| Line 377 \| template<typename T> void ZConstraint<T>::update()
377
378		zconsMols.clear();
379		massOfZConsMols.clear();
380	+	zPos.clear();
381	+	kz.clear();
382
383		unconsMols.clear();
384		massOfUnconsMols.clear();
#	Line 230 \| 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		massOfZConsMols.push_back(molecules[i].getTotalMass());
398
399		molecules[i].getCOM(COM);
#	Line 242 \| Line 406 \| template<typename T> void ZConstraint<T>::update()
406
407		}
408		}
409	+
410	+	//determine the states of z-constraint molecules
411	+	for(int i = 0; i < (int)(zconsMols.size()); i++){
412	+	zconsMols[i]->getCOM(COM);
413	+	if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
414	+	states.push_back(zcsFixed);
415	+	else
416	+	states.push_back(zcsMoving);
417	+	}
418	+
419
420		//The reason to declare fz and indexOfZconsMols as pointer to array is
421		// that we want to make the MPI communication simple
422		if(fz)
423		delete[] fz;
424	+
425	+	if(curZPos)
426	+	delete[] curZPos;
427
428		if(indexOfZConsMols)
429		delete[] indexOfZConsMols;
430
431		if (zconsMols.size() > 0){
432		fz = new double[zconsMols.size()];
433	+	curZPos = new double[zconsMols.size()];
434		indexOfZConsMols = new int[zconsMols.size()];
435
436	<	if(!fz \|\| !indexOfZConsMols){
436	>	if(!fz \|\| !curZPos \|\| !indexOfZConsMols){
437		sprintf( painCave.errMsg,
438		"Memory allocation failure in class Zconstraint\n");
439		painCave.isFatal = 1;
440		simError();
441		}
442
443	<	for(int i = 0; i < zconsMols.size(); i++){
443	>	for(int i = 0; i < (int)(zconsMols.size()); i++){
444		indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
445		}
446
447		}
448		else{
449		fz = NULL;
450	+	curZPos = NULL;
451		indexOfZConsMols = NULL;
452		}
453
454	+	//
455	+	forcePolicy->update();
456	+
457		}
458
459		#endif
460
461	<	/** Function Name: isZConstraintMol
462	<	** Parameter
463	<	** Molecule* mol
464	<	** Return value:
465	<	** -1, if the molecule is not z-constraint molecule,
466	<	** other non-negative values, its index in indexOfAllZConsMols vector
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
468		*/
469
470		template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
#	Line 294 \| Line 477 \| template<typename T> int ZConstraint<T>::isZConstraint
477		index = mol->getGlobalIndex();
478
479		low = 0;
480	<	high = indexOfAllZConsMols.size() - 1;
480	>	high = parameters->size() - 1;
481
482		//Binary Search (we have sorted the array)
483		while(low <= high){
484		mid = (low + high) /2;
485	<	if (indexOfAllZConsMols[mid] == index)
485	>	if ((*parameters)[mid].zconsIndex == index)
486		return mid;
487	<	else if (indexOfAllZConsMols[mid] > index )
487	>	else if ((*parameters)[mid].zconsIndex > index )
488		high = mid -1;
489		else
490		low = mid + 1;
#	Line 310 \| Line 493 \| template<typename T> int ZConstraint<T>::isZConstraint
493		return -1;
494		}
495
313	–	/**
314	–	* Description:
315	–	* Reset the z coordinates
316	–	*/
496		template<typename T> void ZConstraint<T>::integrate(){
497	+
498	+	// creat zconsWriter
499	+	fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);
500
501	+	if(!fzOut){
502	+	sprintf( painCave.errMsg,
503	+	"Memory allocation failure in class Zconstraint\n");
504	+	painCave.isFatal = 1;
505	+	simError();
506	+	}
507	+
508		//zero out the velocities of center of mass of unconstrained molecules
509		//and the velocities of center of mass of every single z-constrained molecueles
510		zeroOutVel();
511	+
512	+	curZconsTime = zconsTime + info->getTime();
513
514		T::integrate();
515
#	Line 330 \| Line 521 \| template<typename T> void ZConstraint<T>::integrate(){
521		*
522		*
523		*
524	<	*/
334	<
335	<
524	>	*/
525		template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
526	+	double zsys;
527	+	double COM[3];
528	+	double force[3];
529	+	double zSysCOMVel;
530
531		T::calcForce(calcPot, calcStress);
532
533	<	if (checkZConsState())
534	<	zeroOutVel();
533	>	if (checkZConsState()){
534	>	zeroOutVel();
535	>	forcePolicy->update();
536	>	}
537	>
538	>	zsys = calcZSys();
539	>	zSysCOMVel = calcSysCOMVel();
540	>	#ifdef IS_MPI
541	>	if(worldRank == 0){
542	>	#endif
543	>	//cout << "---------------------------------------------------------------------" <<endl;
544	>	//cout << "current time: " << info->getTime() << endl;
545	>	//cout << "center of mass at z: " << zsys << endl;
546	>	//cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
547
548	+	#ifdef IS_MPI
549	+	}
550	+	#endif
551	+
552		//do zconstraint force;
553		if (haveFixedZMols())
554		this->doZconstraintForce();
555	<
555	>
556		//use harmonical poteintial to move the molecules to the specified positions
557		if (haveMovingZMols())
558	<	//this->doHarmonic();
559	<
560	<	fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
561	<
558	>	this->doHarmonic();
559	>
560	>	//write out forces and current positions of z-constraint molecules
561	>	if(info->getTime() >= curZconsTime){
562	>	for(int i = 0; i < (int)(zconsMols.size()); i++){
563	>	zconsMols[i]->getCOM(COM);
564	>	curZPos[i] = COM[whichDirection];
565	>
566	>	//if the z-constraint molecule is still moving, just record its force
567	>	if(states[i] == zcsMoving){
568	>	fz[i] = 0;
569	>	Atom** movingZAtoms;
570	>	movingZAtoms = zconsMols[i]->getMyAtoms();
571	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
572	>	movingZAtoms[j]->getFrc(force);
573	>	fz[i] += force[whichDirection];
574	>	}
575	>	}
576	>	}
577	>	fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos);
578	>	curZconsTime += zconsTime;
579	>	}
580	>
581	>	zSysCOMVel = calcSysCOMVel();
582	>	#ifdef IS_MPI
583	>	if(worldRank == 0){
584	>	#endif
585	>	//cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
586	>	#ifdef IS_MPI
587	>	}
588	>	#endif
589		}
590	+
591	+
592	+	/**
593	+	*
594	+	*/
595
596		template<typename T> double ZConstraint<T>::calcZSys()
597		{
#	Line 359 \| Line 600 \| template<typename T> double ZConstraint<T>::calcZSys()
600		double totalMass;
601		double totalMZ_local;
602		double totalMZ;
362	–	double massOfUncons_local;
603		double massOfCurMol;
604		double COM[3];
605
606		totalMass_local = 0;
367	–	totalMass = 0;
607		totalMZ_local = 0;
369	–	totalMZ = 0;
370	–	massOfUncons_local = 0;
371	–
608
609		for(int i = 0; i < nMols; i++){
610		massOfCurMol = molecules[i].getTotalMass();
#	Line 376 \| Line 612 \| template<typename T> double ZConstraint<T>::calcZSys()
612
613		totalMass_local += massOfCurMol;
614		totalMZ_local += massOfCurMol * COM[whichDirection];
615	<
380	<	if(isZConstraintMol(&molecules[i]) == -1){
381	<
382	<	massOfUncons_local += massOfCurMol;
383	<	}
384	<
615	>
616		}
617	+
618
387	–
619		#ifdef IS_MPI
620	<	MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
621	<	MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
622	<	MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
623	<	#else
620	>	MPI_Allreduce(&totalMass_local, &totalMass, 1,
621	>	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
622	>	MPI_Allreduce(&totalMZ_local, &totalMZ, 1,
623	>	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
624	>	#else
625		totalMass = totalMass_local;
626		totalMZ = totalMZ_local;
627	<	totalMassOfUncons = massOfUncons_local;
396	<	#endif
627	>	#endif
628
629		double zsys;
630		zsys = totalMZ / totalMass;
#	Line 412 \| Line 643 \| template<typename T> void ZConstraint<T>::thermalize(
643
644		/**
645		*
415	–	*
416	–	*
646		*/
647
648		template<typename T> void ZConstraint<T>::zeroOutVel(){
#	Line 421 \| Line 650 \| template<typename T> void ZConstraint<T>::zeroOutVel()
650		Atom** fixedZAtoms;
651		double COMvel[3];
652		double vel[3];
653	<
653	>	double zSysCOMVel;
654	>
655		//zero out the velocities of center of mass of fixed z-constrained molecules
656
657	<	for(int i = 0; i < zconsMols.size(); i++){
657	>	for(int i = 0; i < (int)(zconsMols.size()); i++){
658
659	<	if (states[i] == zcsFixed){
659	>	if (states[i] == zcsFixed){
660
661	<	zconsMols[i]->getCOMvel(COMvel);
661	>	zconsMols[i]->getCOMvel(COMvel);
662	>	//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
663	>
664		fixedZAtoms = zconsMols[i]->getMyAtoms();
665	<
665	>
666		for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
667		fixedZAtoms[j]->getVel(vel);
668	<	vel[whichDirection] -= COMvel[whichDirection];
669	<	fixedZAtoms[j]->setVel(vel);
668	>	vel[whichDirection] -= COMvel[whichDirection];
669	>	fixedZAtoms[j]->setVel(vel);
670		}
671	<
671	>
672	>	zconsMols[i]->getCOMvel(COMvel);
673	>	//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
674		}
441	–
442	–	}
675
676	+	}
677	+
678	+	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
679	+
680	+	zSysCOMVel = calcSysCOMVel();
681	+	#ifdef IS_MPI
682	+	if(worldRank == 0){
683	+	#endif
684	+	//cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;
685	+	#ifdef IS_MPI
686	+	}
687	+	#endif
688	+
689		// calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
690		double MVzOfMovingMols_local;
691		double MVzOfMovingMols;
#	Line 450 \| Line 695 \| template<typename T> void ZConstraint<T>::zeroOutVel()
695		MVzOfMovingMols_local = 0;
696		totalMassOfMovingZMols_local = 0;
697
698	<	for(int i =0; i < unconsMols.size(); i++){
698	>	for(int i =0; i < (int)(unconsMols.size()); i++){
699		unconsMols[i]->getCOMvel(COMvel);
700		MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
701		}
702
703	<	for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
459	<
703	>	for(int i = 0; i < (int)(zconsMols.size()); i++){
704		if (states[i] == zcsMoving){
705		zconsMols[i]->getCOMvel(COMvel);
706		MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
707	<	totalMassOfMovingZMols_local += massOfZConsMols[i];
707	>	totalMassOfMovingZMols_local += massOfZConsMols[i];
708		}
709	<
709	>
710		}
711
712		#ifndef IS_MPI
#	Line 478 \| Line 722 \| template<typename T> void ZConstraint<T>::zeroOutVel()
722
723		//modify the velocites of unconstrained molecules
724		Atom** unconsAtoms;
725	<	for(int i = 0; i < unconsMols.size(); i++){
725	>	for(int i = 0; i < (int)(unconsMols.size()); i++){
726
727		unconsAtoms = unconsMols[i]->getMyAtoms();
728		for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
#	Line 491 \| Line 735 \| template<typename T> void ZConstraint<T>::zeroOutVel()
735
736		//modify the velocities of moving z-constrained molecuels
737		Atom** movingZAtoms;
738	<	for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
738	>	for(int i = 0; i < (int)(zconsMols.size()); i++){
739
740		if (states[i] ==zcsMoving){
741
742		movingZAtoms = zconsMols[i]->getMyAtoms();
743	<	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
743	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
744		movingZAtoms[j]->getVel(vel);
745		vel[whichDirection] -= vzOfMovingMols;
746	<	movingZAtoms[j]->setVel(vel);
503	<	}
504	<
746	>	movingZAtoms[j]->setVel(vel);
747		}
748	+
749	+	}
750
751	+	}
752	+
753	+
754	+	zSysCOMVel = calcSysCOMVel();
755	+	#ifdef IS_MPI
756	+	if(worldRank == 0){
757	+	#endif
758	+	//cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;
759	+	#ifdef IS_MPI
760		}
761	+	#endif
762
763		}
764
765	+	/**
766	+	*
767	+	*/
768	+
769		template<typename T> void ZConstraint<T>::doZconstraintForce(){
770
771		Atom** zconsAtoms;
772		double totalFZ;
773		double totalFZ_local;
516	–	double COMvel[3];
774		double COM[3];
775		double force[3];
519	–	double zsys;
776
777	<	int nMovingZMols_local;
522	<	int nMovingZMols;
523	<
524	<	//constrain the molecules which do not reach the specified positions
525	<
526	<	zsys = calcZSys();
527	<	cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;
777	>	//constrain the molecules which do not reach the specified positions
778
779		//Zero Out the force of z-contrained molecules
780		totalFZ_local = 0;
781
782		//calculate the total z-contrained force of fixed z-contrained molecules
783	<	for(int i = 0; i < zconsMols.size(); i++){
784	<
783	>
784	>	//cout << "before zero out z-constraint force on fixed z-constraint molecuels "
785	>	// << "total force is " << calcTotalForce() << endl;
786	>
787	>	for(int i = 0; i < (int)(zconsMols.size()); i++){
788	>
789		if (states[i] == zcsFixed){
790	<
790	>
791		zconsMols[i]->getCOM(COM);
792		zconsAtoms = zconsMols[i]->getMyAtoms();
793
#	Line 541 \| Line 795 \| template<typename T> void ZConstraint<T>::doZconstrain
795		for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
796		zconsAtoms[j]->getFrc(force);
797		fz[i] += force[whichDirection];
798	<	}
798	>	}
799		totalFZ_local += fz[i];
800
801	<	cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
801	>	//cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
802	>	// <<"\tcurrent zpos: " << COM[whichDirection]
803	>	// << "\tcurrent fz: " <<fz[i] << endl;
804
805	+
806		}
807	<
807	>
808		}
809
810	<	//calculate the number of atoms of moving z-constrained molecules
554	<	nMovingZMols_local = 0;
555	<	for(int i = 0; zconsMols.size(); i++){
556	<	if(states[i] == zcsMoving)
557	<	nMovingZMols_local += massOfZConsMols[i];
558	<	}
810	>	//calculate total z-constraint force
811		#ifdef IS_MPI
812		MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
561	–	MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
813		#else
814		totalFZ = totalFZ_local;
564	–	nMovingZMols = nMovingZMols_local;
815		#endif
816
817	+
818	+	// apply negative to fixed z-constrained molecues;
819		force[0]= 0;
820		force[1]= 0;
821		force[2]= 0;
570	–	force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);
822
823	<	//modify the velocites of unconstrained molecules
824	<	for(int i = 0; i < unconsMols.size(); i++){
823	>	for(int i = 0; i < (int)(zconsMols.size()); i++){
824	>
825	>	if (states[i] == zcsFixed){
826	>
827	>	int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
828	>	zconsAtoms = zconsMols[i]->getMyAtoms();
829	>
830	>	for(int j =0; j < nAtomOfCurZConsMol; j++) {
831	>	//force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
832	>	force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
833	>	zconsAtoms[j]->addFrc(force);
834	>	}
835	>
836	>	}
837	>
838	>	}
839	>
840	>	//cout << "after zero out z-constraint force on fixed z-constraint molecuels "
841	>	// << "total force is " << calcTotalForce() << endl;
842	>
843	>
844	>	force[0]= 0;
845	>	force[1]= 0;
846	>	force[2]= 0;
847	>
848	>	//modify the forces of unconstrained molecules
849	>	for(int i = 0; i < (int)(unconsMols.size()); i++){
850
851		Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
852
853	<	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)
853	>	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
854	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
855	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
856		unconsAtoms[j]->addFrc(force);
857	+	}
858
859		}
860
861	<	//modify the velocities of moving z-constrained molecules
862	<	for(int i = 0; i < zconsMols.size(); i++) {
863	<	if (states[i] == zcsMoving){
864	<
865	<	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
861	>	//modify the forces of moving z-constrained molecules
862	>	for(int i = 0; i < (int)(zconsMols.size()); i++) {
863	>	if (states[i] == zcsMoving){
864	>
865	>	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
866
867	<	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)
868	<	movingZAtoms[j]->addFrc(force);
869	<	}
867	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
868	>	//force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
869	>	force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
870	>	movingZAtoms[j]->addFrc(force);
871	>	}
872	>	}
873		}
874	+	// cout << "after substracting z-constraint force from moving molecuels "
875	+	// << "total force is " << calcTotalForce() << endl;
876
877	<	// apply negative to fixed z-constrained molecues;
877	>	}
878	>
879	>	/**
880	>	*
881	>	*
882	>	*/
883	>
884	>	template<typename T> void ZConstraint<T>::doHarmonic(){
885	>	double force[3];
886	>	double harmonicU;
887	>	double harmonicF;
888	>	double COM[3];
889	>	double diff;
890	>	double totalFZ_local;
891	>	double totalFZ;
892	>
893	>	force[0] = 0;
894	>	force[1] = 0;
895	>	force[2] = 0;
896	>
897	>	totalFZ_local = 0;
898	>
899	>	for(int i = 0; i < (int)(zconsMols.size()); i++) {
900	>
901	>	if (states[i] == zcsMoving){
902	>	zconsMols[i]->getCOM(COM);
903	>	// cout << "Moving Molecule\tindex: " << indexOfZConsMols[i]
904	>	// << "\tcurrent zpos: " << COM[whichDirection] << endl;
905	>
906	>	diff = COM[whichDirection] -zPos[i];
907	>
908	>	harmonicU = 0.5 * kz[i] * diff * diff;
909	>	info->lrPot += harmonicU;
910	>
911	>	harmonicF = - kz[i] * diff;
912	>	totalFZ_local += harmonicF;
913	>
914	>	//adjust force
915	>
916	>	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
917	>
918	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
919	>	//force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
920	>	force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
921	>	movingZAtoms[j]->addFrc(force);
922	>	}
923	>	}
924	>
925	>	}
926	>
927	>	#ifndef IS_MPI
928	>	totalFZ = totalFZ_local;
929	>	#else
930	>	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
931	>	#endif
932	>
933	>	//cout << "before substracting harmonic force from moving molecuels "
934	>	// << "total force is " << calcTotalForce() << endl;
935	>
936		force[0]= 0;
937		force[1]= 0;
938		force[2]= 0;
939
940	<	for(int i = 0; i < zconsMols.size(); i++){
940	>	//modify the forces of unconstrained molecules
941	>	for(int i = 0; i < (int)(unconsMols.size()); i++){
942	>
943	>	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
944	>
945	>	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
946	>	//force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
947	>	force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
948	>	unconsAtoms[j]->addFrc(force);
949	>	}
950	>	}
951
952	<	if (states[i] == zcsFixed){
953	<
602	<	int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
603	<	zconsAtoms = zconsMols[i]->getMyAtoms();
604	<
605	<	for(int j =0; j < nAtomOfCurZConsMol; j++) {
606	<	force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
607	<	zconsAtoms[j]->addFrc(force);
608	<	}
609	<
610	<	}
611	<
612	<	}
952	>	//cout << "after substracting harmonic force from moving molecuels "
953	>	// << "total force is " << calcTotalForce() << endl;
954
955		}
956
957	+	/**
958	+	*
959	+	*/
960	+
961		template<typename T> bool ZConstraint<T>::checkZConsState(){
962		double COM[3];
963		double diff;
964
965	<	bool changed;
965	>	int changed_local;
966	>	int changed;
967
968	<	changed = false;
968	>	changed_local = 0;
969
970	<	for(int i =0; i < zconsMols.size(); i++){
970	>	for(int i =0; i < (int)(zconsMols.size()); i++){
971
972		zconsMols[i]->getCOM(COM);
973	<	diff = fabs(COM[whichDirection] - ZPos[i]);
974	<	if ( diff <= ztol && states[i] == zcsMoving){
973	>	diff = fabs(COM[whichDirection] - zPos[i]);
974	>	if ( diff <= zconsTol && states[i] == zcsMoving){
975		states[i] = zcsFixed;
976	<	changed = true;
976	>	changed_local = 1;
977		}
978	<	else if ( diff > ztol && states[i] == zcsFixed){
978	>	else if ( diff > zconsTol && states[i] == zcsFixed){
979		states[i] = zcsMoving;
980	<	changed = true;
980	>	changed_local = 1;
981		}
982
983		}
984
985	<	return changed;
985	>	#ifndef IS_MPI
986	>	changed =changed_local;
987	>	#else
988	>	MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
989	>	#endif
990	>
991	>	return (changed > 0);
992	>
993		}
994
995		template<typename T> bool ZConstraint<T>::haveFixedZMols(){
643	–	for(int i = 0; i < zconsMols.size(); i++)
644	–	if (states[i] == zcsFixed)
645	–	return true;
996
997	<	return false;
997	>	int havingFixed_local;
998	>	int havingFixed;
999	>
1000	>	havingFixed_local = 0;
1001	>
1002	>	for(int i = 0; i < (int)(zconsMols.size()); i++)
1003	>	if (states[i] == zcsFixed){
1004	>	havingFixed_local = 1;
1005	>	break;
1006	>	}
1007	>
1008	>	#ifndef IS_MPI
1009	>	havingFixed = havingFixed_local;
1010	>	#else
1011	>	MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1012	>	#endif
1013	>
1014	>	return (havingFixed > 0);
1015		}
1016
1017
#	Line 652 \| Line 1019 \| template<typename T> bool ZConstraint<T>::haveMovingZM
1019		*
1020		*/
1021		template<typename T> bool ZConstraint<T>::haveMovingZMols(){
655	–	for(int i = 0; i < zconsMols.size(); i++)
656	–	if (states[i] == zcsMoving)
657	–	return true;
1022
1023	<	return false;
1023	>	int havingMoving_local;
1024	>	int havingMoving;
1025	>
1026	>	havingMoving_local = 0;
1027	>
1028	>	for(int i = 0; i < (int)(zconsMols.size()); i++)
1029	>	if (states[i] == zcsMoving){
1030	>	havingMoving_local = 1;
1031	>	break;
1032	>	}
1033	>
1034	>	#ifndef IS_MPI
1035	>	havingMoving = havingMoving_local;
1036	>	#else
1037	>	MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
1038	>	#endif
1039	>
1040	>	return (havingMoving > 0);
1041
1042	<	}
1042	>	}
1043	>
1044	>	/**
1045	>	*
1046	>	*/
1047	>
1048	>	template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
1049	>	{
1050	>	double MVzOfMovingMols_local;
1051	>	double MVzOfMovingMols;
1052	>	double totalMassOfMovingZMols_local;
1053	>	double totalMassOfMovingZMols;
1054	>	double COMvel[3];
1055	>
1056	>	MVzOfMovingMols_local = 0;
1057	>	totalMassOfMovingZMols_local = 0;
1058	>
1059	>	for(int i =0; i < unconsMols.size(); i++){
1060	>	unconsMols[i]->getCOMvel(COMvel);
1061	>	MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
1062	>	}
1063	>
1064	>	for(int i = 0; i < zconsMols.size(); i++){
1065	>
1066	>	if (states[i] == zcsMoving){
1067	>	zconsMols[i]->getCOMvel(COMvel);
1068	>	MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
1069	>	totalMassOfMovingZMols_local += massOfZConsMols[i];
1070	>	}
1071	>
1072	>	}
1073	>
1074	>	#ifndef IS_MPI
1075	>	MVzOfMovingMols = MVzOfMovingMols_local;
1076	>	totalMassOfMovingZMols = totalMassOfMovingZMols_local;
1077	>	#else
1078	>	MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1079	>	MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1080	>	#endif
1081	>
1082	>	double vzOfMovingMols;
1083	>	vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
1084	>
1085	>	return vzOfMovingMols;
1086	>	}
1087	>
1088	>	/**
1089	>	*
1090	>	*/
1091	>
1092	>	template<typename T> double ZConstraint<T>::calcSysCOMVel()
1093	>	{
1094	>	double COMvel[3];
1095	>	double tempMVz_local;
1096	>	double tempMVz;
1097	>	double massOfZCons_local;
1098	>	double massOfZCons;
1099	>
1100	>
1101	>	tempMVz_local = 0;
1102	>
1103	>	for(int i =0 ; i < nMols; i++){
1104	>	molecules[i].getCOMvel(COMvel);
1105	>	tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
1106	>	}
1107	>
1108	>	massOfZCons_local = 0;
1109	>
1110	>	for(int i = 0; i < (int)(massOfZConsMols.size()); i++){
1111	>	massOfZCons_local += massOfZConsMols[i];
1112	>	}
1113	>	#ifndef IS_MPI
1114	>	massOfZCons = massOfZCons_local;
1115	>	tempMVz = tempMVz_local;
1116	>	#else
1117	>	MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1118	>	MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1119	>	#endif
1120	>
1121	>	return tempMVz /(totalMassOfUncons + massOfZCons);
1122	>	}
1123	>
1124	>	/**
1125	>	*
1126	>	*/
1127	>
1128	>	template<typename T> double ZConstraint<T>::calcTotalForce(){
1129	>
1130	>	double force[3];
1131	>	double totalForce_local;
1132	>	double totalForce;
1133	>
1134	>	totalForce_local = 0;
1135	>
1136	>	for(int i = 0; i < nAtoms; i++){
1137	>	atoms[i]->getFrc(force);
1138	>	totalForce_local += force[whichDirection];
1139	>	}
1140	>
1141	>	#ifndef IS_MPI
1142	>	totalForce = totalForce_local;
1143	>	#else
1144	>	MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1145	>	#endif
1146	>
1147	>	return totalForce;
1148	>
1149	>	}
1150	>
1151	>	/**
1152	>	*
1153	>	*/
1154	>
1155	>	template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
1156	>	//calculate the number of atoms of moving z-constrained molecules
1157	>	int nMovingZAtoms_local;
1158	>	int nMovingZAtoms;
1159	>
1160	>	nMovingZAtoms_local = 0;
1161	>	for(int i = 0; i < (int)((zconsIntegrator->zconsMols).size()); i++)
1162	>	if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1163	>	nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1164	>
1165	>	#ifdef IS_MPI
1166	>	MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1167	>	#else
1168	>	nMovingZAtoms = nMovingZAtoms_local;
1169	>	#endif
1170	>	totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1171	>	}
1172	>
1173	>	template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1174	>	return totalForce / mol->getNAtoms();
1175	>	}
1176	>
1177	>	template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, double totalForce){
1178	>	return totalForce / totNumOfMovingAtoms;
1179	>	}
1180	>
1181	>	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1182	>	return totalForce / mol->getNAtoms();
1183	>	}
1184	>
1185	>	template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, double totalForce){
1186	>	return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
1187	>	}
1188	>
1189	>	/**
1190	>	*
1191	>	*/
1192	>
1193	>	template<typename T> void ZConstraint<T>::PolicyByMass::update(){
1194	>	//calculate the number of atoms of moving z-constrained molecules
1195	>	double massOfMovingZAtoms_local;
1196	>	double massOfMovingZAtoms;
1197	>
1198	>	massOfMovingZAtoms_local = 0;
1199	>	for(int i = 0; i < (int)((zconsIntegrator->zconsMols).size()); i++)
1200	>	if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1201	>	massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1202	>
1203	>	#ifdef IS_MPI
1204	>	MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1205	>	#else
1206	>	massOfMovingZAtoms = massOfMovingZAtoms_local;
1207	>	#endif
1208	>	totMassOfMovingAtoms = massOfMovingZAtoms + zconsIntegrator->totalMassOfUncons;
1209	>	}
1210	>
1211	>	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1212	>	return totalForce * atom->getMass() / mol->getTotalMass();
1213	>	}
1214	>
1215	>	template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols( Atom* atom, double totalForce){
1216	>	return totalForce * atom->getMass() / totMassOfMovingAtoms;
1217	>	}
1218	>
1219	>	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1220	>	return totalForce * atom->getMass() / mol->getTotalMass();
1221	>	}
1222	>
1223	>	template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, double totalForce){
1224	>	return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
1225	>	}
1226	>

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents): Revision 676 by tim, Mon Aug 11 19:40:06 2003 UTC vs. Revision 830 by gezelter, Tue Oct 28 16:20:28 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 676 by tim, Mon Aug 11 19:40:06 2003 UTC vs.
Revision 830 by gezelter, Tue Oct 28 16:20:28 2003 UTC