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

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 699 by tim, Fri Aug 15 19:24:13 2003 UTC vs.
Revision 702 by tim, Wed Aug 20 14:50:32 2003 UTC

#	Line 3 \| Line 3 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
3		#include <cmath>
4		template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff)
5		: T(theInfo, the_ff), fz(NULL), curZPos(NULL),
6	<	indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0)
6	>	indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0)
7		{
8
9		//get properties from SimInfo
#	Line 40 \| Line 40 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
40		}
41		else{
42		policy = dynamic_cast<StringData*>(data);
43	<
44	<	if(!policy){
43	>
44	>	if(!policy){
45		sprintf( painCave.errMsg,
46		"ZConstraint Error: Convertion from GenericData to StringData failure, "
47		"average force substraction policy is used\n");
#	Line 49 \| Line 49 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
49		simError();
50
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{
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 Warning: unknown force substraction policy, "
61		"average force substraction policy is used\n");
62		painCave.isFatal = 0;
63		simError();
64	<	}
65	<	}
64	>	}
65	>	}
66		}
67	–
67
68	<
68	>
69		//retrieve sample time of z-contraint
70		data = info->getProperty(ZCONSTIME_ID);
71
#	Line 110 \| Line 109 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
109		}
110		else{
111
112	<	filename = dynamic_cast<StringData*>(data);
112	>	filename = dynamic_cast<StringData*>(data);
113
114		if(!filename){
115
#	Line 124 \| Line 123 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
123		this->zconsOutput = filename->getData();
124		}
125
127	–
126		}
127
128		//retrieve tolerance for z-constraint molecuels
#	Line 154 \| Line 152 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
152		}
153
154		}
155	<
155	>
156		//retrieve index of z-constraint molecules
157		data = info->getProperty(ZCONSPARADATA_ID);
158		if(!data) {
#	Line 195 \| Line 193 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
193		"ZConstraint error: index is out of range\n");
194		painCave.isFatal = 1;
195		simError();
196	<	}
196	>	}
197
198		maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
199
#	Line 216 \| Line 214 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
214		//its initial z coordinate will be used as default
215		for(int i = 0; i < parameters->size(); i++){
216
217	<	if(!(*parameters)[i].havingZPos){
220	<
217	>	if(!(*parameters)[i].havingZPos){
218		#ifndef IS_MPI
219	<	for(int j = 0; j < nMols; j++){
220	<	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
221	<	molecules[j].getCOM(COM);
222	<	break;
226	<	}
219	>	for(int j = 0; j < nMols; j++){
220	>	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
221	>	molecules[j].getCOM(COM);
222	>	break;
223		}
224	+	}
225		#else
226		//query which processor current zconstraint molecule belongs to
227	<	int *MolToProcMap;
228	<	int whichNode;
229	<	double initZPos;
230	<	MolToProcMap = mpiSim->getMolToProcMap();
231	<	whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
232	<
233	<	//broadcast the zpos of current z-contraint molecule
234	<	//the node which contain this
227	>	int *MolToProcMap;
228	>	int whichNode;
229	>	double initZPos;
230	>	MolToProcMap = mpiSim->getMolToProcMap();
231	>	whichNode = MolToProcMap[(*parameters)[i].zconsIndex];
232	>
233	>	//broadcast the zpos of current z-contraint molecule
234	>	//the node which contain this
235
236	<	if (worldRank == whichNode ){
237	<
238	<	for(int j = 0; j < nMols; j++)
239	<	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
240	<	molecules[j].getCOM(COM);
241	<	break;
242	<	}
243	<
244	<	}
236	>	if (worldRank == whichNode ){
237	>
238	>	for(int j = 0; j < nMols; j++)
239	>	if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){
240	>	molecules[j].getCOM(COM);
241	>	break;
242	>	}
243	>
244	>	}
245
246	<	MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);
246	>	MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);
247		#endif
248
249	<	(*parameters)[i].zPos = COM[whichDirection];
249	>	(*parameters)[i].zPos = COM[whichDirection];
250
251	<	sprintf( painCave.errMsg,
251	>	sprintf( painCave.errMsg,
252		"ZConstraint warningr: Does not specify zpos for z-constraint molecule "
253		"initial z coornidate will be used \n");
254	<	painCave.isFatal = 0;
255	<	simError();
256	<
257	<	}
258	<	}
259	<
254	>	painCave.isFatal = 0;
255	>	simError();
256	>
257	>	}
258	>	}
259	>
260		}//end if (!zConsParaData)
261		}//end if (!data)
262
#	Line 279 \| Line 276 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
276		massOfZConsMols.push_back(molecules[i].getTotalMass());
277
278		zPos.push_back((*parameters)[searchResult].zPos);
279	<	cout << "index: "<< (parameters)[searchResult].zconsIndex <<"\tzPos = " << (parameters)[searchResult].zPos << endl;
280	<	kz.push_back((parameters)[searchResult]. kRatio zForceConst);
279	>	cout << "index: "<< (*parameters)[searchResult].zconsIndex
280	>	<<"\tzPos = " << (*parameters)[searchResult].zPos << endl;
281	>	kz.push_back((parameters)[searchResult]. kRatio zForceConst);
282
283		molecules[i].getCOM(COM);
284		}
#	Line 308 \| Line 306 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
306		for(int i = 0; i < zconsMols.size(); i++){
307		indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
308
309	<	zconsMols[i]->getCOM(COM);
309	>	zconsMols[i]->getCOM(COM);
310		if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
311	<	states.push_back(zcsFixed);
312	<	else
313	<	states.push_back(zcsMoving);
311	>	states.push_back(zcsFixed);
312	>	else
313	>	states.push_back(zcsMoving);
314		}
315
316		#endif
#	Line 327 \| Line 325 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
325		#ifndef IS_MPI
326		totalMassOfUncons = totalMassOfUncons_local;
327		#else
328	<	MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
328	>	MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1,
329	>	MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
330		#endif
331
332
#	Line 340 \| Line 339 \| template<typename T> ZConstraint<T>::ZConstraint(SimIn
339		#ifndef IS_MPI
340		totNumOfUnconsAtoms = nUnconsAtoms_local;
341		#else
342	<	MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
342	>	MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1,
343	>	MPI_INT,MPI_SUM, MPI_COMM_WORLD);
344		#endif
345
346		// creat zconsWriter
#	Line 369 \| Line 369 \| template<typename T> ZConstraint<T>::~ZConstraint()
369
370		if(fzOut)
371		delete fzOut;
372	<
372	>
373		if(forcePolicy)
374		delete forcePolicy;
375		}
376
377	+
378	+	/**
379	+	*
380	+	*/
381	+
382		#ifdef IS_MPI
383		template<typename T> void ZConstraint<T>::update()
384		{
#	Line 398 \| Line 403 \| template<typename T> void ZConstraint<T>::update()
403
404		zconsMols.push_back(&molecules[i]);
405		zPos.push_back((*parameters)[index].zPos);
406	<	kz.push_back((parameters)[index].kRatio zForceConst);
407	<
406	>	kz.push_back((parameters)[index].kRatio zForceConst);
407	>
408		massOfZConsMols.push_back(molecules[i].getTotalMass());
409
410		molecules[i].getCOM(COM);
#	Line 415 \| Line 420 \| template<typename T> void ZConstraint<T>::update()
420
421		//determine the states of z-constraint molecules
422		for(int i = 0; i < zconsMols.size(); i++){
423	<	zconsMols[i]->getCOM(COM);
423	>	zconsMols[i]->getCOM(COM);
424		if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
425	<	states.push_back(zcsFixed);
426	<	else
427	<	states.push_back(zcsMoving);
425	>	states.push_back(zcsFixed);
426	>	else
427	>	states.push_back(zcsMoving);
428		}
429
430
#	Line 427 \| Line 432 \| template<typename T> void ZConstraint<T>::update()
432		// that we want to make the MPI communication simple
433		if(fz)
434		delete[] fz;
435	<
435	>
436		if(curZPos)
437		delete[] curZPos;
438
#	Line 436 \| Line 441 \| template<typename T> void ZConstraint<T>::update()
441
442		if (zconsMols.size() > 0){
443		fz = new double[zconsMols.size()];
444	<	curZPos = new double[zconsMols.size()];
444	>	curZPos = new double[zconsMols.size()];
445		indexOfZConsMols = new int[zconsMols.size()];
446
447		if(!fz \|\| !curZPos \|\| !indexOfZConsMols){
#	Line 453 \| Line 458 \| template<typename T> void ZConstraint<T>::update()
458		}
459		else{
460		fz = NULL;
461	<	curZPos = NULL;
461	>	curZPos = NULL;
462		indexOfZConsMols = NULL;
463		}
464	<
464	>
465		//
466		forcePolicy->update();
467
#	Line 464 \| Line 469 \| template<typename T> void ZConstraint<T>::update()
469
470		#endif
471
472	<	/** Function Name: isZConstraintMol
473	<	** Parameter
474	<	** Molecule* mol
475	<	** Return value:
476	<	** -1, if the molecule is not z-constraint molecule,
477	<	** other non-negative values, its index in indexOfAllZConsMols vector
472	>	/**
473	>	* Function Name: isZConstraintMol
474	>	* Parameter
475	>	* Molecule* mol
476	>	* Return value:
477	>	* -1, if the molecule is not z-constraint molecule,
478	>	* other non-negative values, its index in indexOfAllZConsMols vector
479		*/
480
481		template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
#	Line 519 \| Line 525 \| template<typename T> void ZConstraint<T>::calcForce(in
525		double zsys;
526		double COM[3];
527		double force[3];
528	+	double zSysCOMVel;
529
530		T::calcForce(calcPot, calcStress);
531
532		if (checkZConsState()){
533		zeroOutVel();
534	<	forcePolicy->update();
534	>	forcePolicy->update();
535		}
536	+
537		zsys = calcZSys();
538	<	cout << "---------------------------------------------------------------------" <<endl;
539	<	cout << "current time: " << info->getTime() << endl;
540	<	cout << "center of mass at z: " << zsys << endl;
541	<	//cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
542	<	cout << "before calcForce, the COMVel of system is " << calcSysCOMVel() <<endl;
538	>	zSysCOMVel = calcSysCOMVel();
539	>	#ifdef IS_MPI
540	>	if(worldRank == 0){
541	>	#endif
542	>	cout << "---------------------------------------------------------------------" <<endl;
543	>	cout << "current time: " << info->getTime() << endl;
544	>	cout << "center of mass at z: " << zsys << endl;
545	>	cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
546
547	<	//cout << "before doZConstraintForce, totalForce is " << calcTotalForce() << endl;
547	>	#ifdef IS_MPI
548	>	}
549	>	#endif
550
551		//do zconstraint force;
552		if (haveFixedZMols())
#	Line 543 \| Line 556 \| template<typename T> void ZConstraint<T>::calcForce(in
556		if (haveMovingZMols())
557		this->doHarmonic();
558
546	–	//cout << "after doHarmonic, totalForce is " << calcTotalForce() << endl;
547	–
559		//write out forces and current positions of z-constraint molecules
560	<	if(info->getTime() >= curZconsTime){
561	<	for(int i = 0; i < zconsMols.size(); i++){
560	>	if(info->getTime() >= curZconsTime){
561	>	for(int i = 0; i < zconsMols.size(); i++){
562		zconsMols[i]->getCOM(COM);
563	<	curZPos[i] = COM[whichDirection];
563	>	curZPos[i] = COM[whichDirection];
564
565	<	//if the z-constraint molecule is still moving, just record its force
566	<	if(states[i] == zcsMoving){
565	>	//if the z-constraint molecule is still moving, just record its force
566	>	if(states[i] == zcsMoving){
567		fz[i] = 0;
568	<	Atom** movingZAtoms;
569	<	movingZAtoms = zconsMols[i]->getMyAtoms();
570	<	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
568	>	Atom** movingZAtoms;
569	>	movingZAtoms = zconsMols[i]->getMyAtoms();
570	>	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
571		movingZAtoms[j]->getFrc(force);
572		fz[i] += force[whichDirection];
573	<	}
574	<	}
575	<	}
573	>	}
574	>	}
575	>	}
576		fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos);
577	<	curZconsTime += zconsTime;
577	>	curZconsTime += zconsTime;
578		}
579	<
580	<	//cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
581	<	cout << "after calcForce, the COMVel of system is " << calcSysCOMVel() <<endl;
579	>
580	>	zSysCOMVel = calcSysCOMVel();
581	>	#ifdef IS_MPI
582	>	if(worldRank == 0){
583	>	#endif
584	>	cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
585	>	#ifdef IS_MPI
586	>	}
587	>	#endif
588	>
589		}
590	+
591	+
592	+	/**
593	+	*
594	+	*/
595
596		template<typename T> double ZConstraint<T>::calcZSys()
597		{
#	Line 579 \| Line 602 \| template<typename T> double ZConstraint<T>::calcZSys()
602		double totalMZ;
603		double massOfCurMol;
604		double COM[3];
605	<
605	>
606		totalMass_local = 0;
607		totalMZ_local = 0;
608
#	Line 594 \| Line 617 \| template<typename T> double ZConstraint<T>::calcZSys()
617
618
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);
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;
#	Line 618 \| Line 643 \| template<typename T> void ZConstraint<T>::thermalize(
643
644		/**
645		*
621	–	*
622	–	*
646		*/
647
648		template<typename T> void ZConstraint<T>::zeroOutVel(){
#	Line 627 \| Line 650 \| template<typename T> void ZConstraint<T>::zeroOutVel()
650		Atom** fixedZAtoms;
651		double COMvel[3];
652		double vel[3];
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++){
658
659	<	if (states[i] == zcsFixed){
659	>	if (states[i] == zcsFixed){
660
661	<	zconsMols[i]->getCOMvel(COMvel);
662	<	//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
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
672	<	zconsMols[i]->getCOMvel(COMvel);
673	<	//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
672	>	zconsMols[i]->getCOMvel(COMvel);
673	>	//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
674		}
675	<
675	>
676		}
677
678	<	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
678	>	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
679
680	+	zSysCOMVel = calcSysCOMVel();
681		#ifdef IS_MPI
682	<	if (worldRank == 0){
682	>	if(worldRank == 0){
683		#endif
684	<	cout << "before resetting the COMVel of sytem is " << calcSysCOMVel() << endl;
684	>	cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;
685		#ifdef IS_MPI
686		}
687		#endif
688	<
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 679 \| Line 704 \| template<typename T> void ZConstraint<T>::zeroOutVel()
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 715 \| Line 740 \| template<typename T> void ZConstraint<T>::zeroOutVel()
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);
747	<	}
748	<
746	>	movingZAtoms[j]->setVel(vel);
747	>	}
748	>
749		}
750
751		}
752
753	+
754	+	zSysCOMVel = calcSysCOMVel();
755		#ifdef IS_MPI
756	<	if (worldRank == 0){
756	>	if(worldRank == 0){
757		#endif
758	<	cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() <<endl;
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;
#	Line 746 \| Line 777 \| template<typename T> void ZConstraint<T>::doZconstrain
777
778
779
780	<	//constrain the molecules which do not reach the specified positions
780	>	//constrain the molecules which do not reach the specified positions
781
782		//Zero Out the force of z-contrained molecules
783		totalFZ_local = 0;
784
785		//calculate the total z-contrained force of fixed z-contrained molecules
786
756	–	#ifdef IS_MPI
757	–	if (worldRank == 0){
758	–	#endif
759	–	cout << "Fixed Molecules" << endl;
760	–	#ifdef IS_MPI
761	–	}
762	–	#endif
763	–
787		for(int i = 0; i < zconsMols.size(); i++){
788	<
788	>
789		if (states[i] == zcsFixed){
790	<
790	>
791		zconsMols[i]->getCOM(COM);
792		zconsAtoms = zconsMols[i]->getMyAtoms();
793
#	Line 772 \| 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]
802	<	<<"\tcurrent zpos: " << COM[whichDirection]
803	<	<< "\tcurrent fz: " <<fz[i] << endl;
801	>	cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i]
802	>	<<"\tcurrent zpos: " << COM[whichDirection]
803	>	<< "\tcurrent fz: " <<fz[i] << endl;
804
805		}
806	<
806	>
807		}
808
809		//calculate total z-constraint force
#	Line 790 \| Line 813 \| template<typename T> void ZConstraint<T>::doZconstrain
813		totalFZ = totalFZ_local;
814		#endif
815
816	<
816	>
817		// apply negative to fixed z-constrained molecues;
818		force[0]= 0;
819		force[1]= 0;
#	Line 799 \| Line 822 \| template<typename T> void ZConstraint<T>::doZconstrain
822		for(int i = 0; i < zconsMols.size(); i++){
823
824		if (states[i] == zcsFixed){
825	<
825	>
826		int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
827		zconsAtoms = zconsMols[i]->getMyAtoms();
828
829		for(int j =0; j < nAtomOfCurZConsMol; j++) {
830	<	force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
830	>	force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
831		//force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
832		zconsAtoms[j]->addFrc(force);
833		}
834	<
834	>
835		}
836	<
836	>
837		}
838
839		//cout << "after zero out z-constraint force on fixed z-constraint molecuels "
840	<	// << "total force is " << calcTotalForce() << endl;
840	>	// << "total force is " << calcTotalForce() << endl;
841
842		//calculate the number of atoms of moving z-constrained molecules
843		int nMovingZAtoms_local;
844		int nMovingZAtoms;
845	<
845	>
846		nMovingZAtoms_local = 0;
847		for(int i = 0; i < zconsMols.size(); i++)
848		if(states[i] == zcsMoving)
849	<	nMovingZAtoms_local += zconsMols[i]->getNAtoms();
849	>	nMovingZAtoms_local += zconsMols[i]->getNAtoms();
850
851		#ifdef IS_MPI
852	<	MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
852	>	MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1,
853	>	MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
854		#else
855		nMovingZAtoms = nMovingZAtoms_local;
856		#endif
#	Line 851 \| Line 875 \| template<typename T> void ZConstraint<T>::doZconstrain
875		//modify the forces of moving z-constrained molecules
876		for(int i = 0; i < zconsMols.size(); i++) {
877		if (states[i] == zcsMoving){
878	<
878	>
879		Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
880
881		for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
#	Line 863 \| Line 887 \| template<typename T> void ZConstraint<T>::doZconstrain
887		}
888
889		//cout << "after substracting z-constraint force from moving molecuels "
890	<	// << "total force is " << calcTotalForce() << endl;
890	>	// << "total force is " << calcTotalForce() << endl;
891
892		}
893
894	+	/**
895	+	*
896	+	*
897	+	*/
898	+
899	+	template<typename T> void ZConstraint<T>::doHarmonic(){
900	+	double force[3];
901	+	double harmonicU;
902	+	double harmonicF;
903	+	double COM[3];
904	+	double diff;
905	+	double totalFZ_local;
906	+	double totalFZ;
907	+
908	+	force[0] = 0;
909	+	force[1] = 0;
910	+	force[2] = 0;
911	+
912	+	totalFZ_local = 0;
913	+
914	+	for(int i = 0; i < zconsMols.size(); i++) {
915	+
916	+	if (states[i] == zcsMoving){
917	+	zconsMols[i]->getCOM(COM);
918	+	cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
919	+
920	+	diff = COM[whichDirection] -zPos[i];
921	+
922	+	harmonicU = 0.5 * kz[i] * diff * diff;
923	+	info->lrPot += harmonicU;
924	+
925	+	harmonicF = - kz[i] * diff;
926	+	totalFZ_local += harmonicF;
927	+
928	+	//adjust force
929	+
930	+	Atom** movingZAtoms = zconsMols[i]->getMyAtoms();
931	+
932	+	for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
933	+	force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
934	+	//force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
935	+	movingZAtoms[j]->addFrc(force);
936	+	}
937	+	}
938	+
939	+	}
940	+
941	+	#ifndef IS_MPI
942	+	totalFZ = totalFZ_local;
943	+	#else
944	+	MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
945	+	#endif
946	+
947	+	force[0]= 0;
948	+	force[1]= 0;
949	+	force[2]= 0;
950	+
951	+	//modify the forces of unconstrained molecules
952	+	for(int i = 0; i < unconsMols.size(); i++){
953	+
954	+	Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
955	+
956	+	for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){
957	+	force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
958	+	//force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
959	+	unconsAtoms[j]->addFrc(force);
960	+	}
961	+	}
962	+
963	+	}
964	+
965	+	/**
966	+	*
967	+	*/
968	+
969		template<typename T> bool ZConstraint<T>::checkZConsState(){
970		double COM[3];
971		double diff;
972
973		int changed_local;
974		int changed;
975	<
975	>
976		changed_local = 0;
977
978		for(int i =0; i < zconsMols.size(); i++){
#	Line 882 \| Line 981 \| template<typename T> bool ZConstraint<T>::checkZConsSt
981		diff = fabs(COM[whichDirection] - zPos[i]);
982		if ( diff <= zconsTol && states[i] == zcsMoving){
983		states[i] = zcsFixed;
984	<	changed_local = 1;
984	>	changed_local = 1;
985		}
986		else if ( diff > zconsTol && states[i] == zcsFixed){
987		states[i] = zcsMoving;
988	<	changed_local = 1;
988	>	changed_local = 1;
989		}
990
991		}
#	Line 910 \| Line 1009 \| template<typename T> bool ZConstraint<T>::haveFixedZMo
1009		for(int i = 0; i < zconsMols.size(); i++)
1010		if (states[i] == zcsFixed){
1011		havingFixed_local = 1;
1012	<	break;
1012	>	break;
1013		}
1014
1015		#ifndef IS_MPI
#	Line 936 \| Line 1035 \| template<typename T> bool ZConstraint<T>::haveMovingZM
1035		for(int i = 0; i < zconsMols.size(); i++)
1036		if (states[i] == zcsMoving){
1037		havingMoving_local = 1;
1038	<	break;
1038	>	break;
1039		}
1040
1041		#ifndef IS_MPI
#	Line 950 \| Line 1049 \| template<typename T> bool ZConstraint<T>::haveMovingZM
1049		}
1050
1051		/**
1052	<	*
1053	<	*
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();
1052	>	*
1053	>	*/
1054
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	–
1055		template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
1056		{
1057		double MVzOfMovingMols_local;
#	Line 1050 \| Line 1073 \| template<typename T> double ZConstraint<T>::calcMoving
1073		if (states[i] == zcsMoving){
1074		zconsMols[i]->getCOMvel(COMvel);
1075		MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];
1076	<	totalMassOfMovingZMols_local += massOfZConsMols[i];
1076	>	totalMassOfMovingZMols_local += massOfZConsMols[i];
1077		}
1078	<
1078	>
1079		}
1080
1081		#ifndef IS_MPI
#	Line 1069 \| Line 1092 \| template<typename T> double ZConstraint<T>::calcMoving
1092		return vzOfMovingMols;
1093		}
1094
1095	+	/**
1096	+	*
1097	+	*/
1098
1099		template<typename T> double ZConstraint<T>::calcSysCOMVel()
1100		{
#	Line 1083 \| Line 1109 \| template<typename T> double ZConstraint<T>::calcSysCOM
1109
1110		for(int i =0 ; i < nMols; i++){
1111		molecules[i].getCOMvel(COMvel);
1112	<	tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
1112	>	tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
1113		}
1114
1115		massOfZCons_local = 0;
1116	<
1116	>
1117		for(int i = 0; i < massOfZConsMols.size(); i++){
1118		massOfZCons_local += massOfZConsMols[i];
1119		}
#	Line 1102 \| Line 1128 \| template<typename T> double ZConstraint<T>::calcSysCOM
1128		return tempMVz /(totalMassOfUncons + massOfZCons);
1129		}
1130
1131	+	/**
1132	+	*
1133	+	*/
1134	+
1135		template<typename T> double ZConstraint<T>::calcTotalForce(){
1136
1137		double force[3];
#	Line 1133 \| Line 1163 \| template<typename T> void ZConstraint<T>::PolicyByNumb
1163		//calculate the number of atoms of moving z-constrained molecules
1164		int nMovingZAtoms_local;
1165		int nMovingZAtoms;
1166	<
1166	>
1167		nMovingZAtoms_local = 0;
1168		for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1169		if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1170	<	nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1170	>	nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
1171
1172		#ifdef IS_MPI
1173		MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#	Line 1147 \| Line 1177 \| template<typename T> void ZConstraint<T>::PolicyByNumb
1177		totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
1178		}
1179
1180	<	template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1180	>	template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
1181		return totalForce / mol->getNAtoms();
1182		}
1183
#	Line 1171 \| Line 1201 \| template<typename T> void ZConstraint<T>::PolicyByMass
1201		//calculate the number of atoms of moving z-constrained molecules
1202		double massOfMovingZAtoms_local;
1203		double massOfMovingZAtoms;
1204	<
1204	>
1205		massOfMovingZAtoms_local = 0;
1206		for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
1207		if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
1208	<	massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1208	>	massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
1209
1210		#ifdef IS_MPI
1211		MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents): Revision 699 by tim, Fri Aug 15 19:24:13 2003 UTC vs. Revision 702 by tim, Wed Aug 20 14:50:32 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 699 by tim, Fri Aug 15 19:24:13 2003 UTC vs.
Revision 702 by tim, Wed Aug 20 14:50:32 2003 UTC