[ViewVC] Diff of: OpenMD/branches/development/src/brains/SimInfo.cpp

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 945 by gezelter, Tue Apr 25 02:09:01 2006 UTC vs.
Revision 1277 by gezelter, Mon Jul 14 12:35:58 2008 UTC

#	Line 53 \| Line 53
53		#include "brains/SimInfo.hpp"
54		#include "math/Vector3.hpp"
55		#include "primitives/Molecule.hpp"
56	+	#include "primitives/StuntDouble.hpp"
57		#include "UseTheForce/fCutoffPolicy.h"
58		#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
59		#include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
60		#include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
61		#include "UseTheForce/doForces_interface.h"
62	+	#include "UseTheForce/DarkSide/neighborLists_interface.h"
63		#include "UseTheForce/DarkSide/electrostatic_interface.h"
64		#include "UseTheForce/DarkSide/switcheroo_interface.h"
65		#include "utils/MemoryUtils.hpp"
#	Line 66 \| Line 68
68		#include "io/ForceFieldOptions.hpp"
69		#include "UseTheForce/ForceField.hpp"
70
71	+
72		#ifdef IS_MPI
73		#include "UseTheForce/mpiComponentPlan.h"
74		#include "UseTheForce/DarkSide/simParallel_interface.h"
#	Line 87 \| Line 90 \| namespace oopse {
90		ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
91		nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
92		nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
93	<	nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0),
94	<	nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0),
95	<	sman_(NULL), fortranInitialized_(false) {
93	>	nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0),
94	>	nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0),
95	>	nConstraints_(0), sman_(NULL), fortranInitialized_(false),
96	>	calcBoxDipole_(false), useAtomicVirial_(true) {
97
98	+
99		MoleculeStamp* molStamp;
100		int nMolWithSameStamp;
101		int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
#	Line 98 \| Line 103 \| namespace oopse {
103		CutoffGroupStamp* cgStamp;
104		RigidBodyStamp* rbStamp;
105		int nRigidAtoms = 0;
106	+
107		std::vector<Component*> components = simParams->getComponents();
108
109		for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
#	Line 152 \| Line 158 \| namespace oopse {
158		+ nGlobalRigidBodies_;
159
160		nGlobalMols_ = molStampIds_.size();
155	–
156	–	#ifdef IS_MPI
161		molToProcMap_.resize(nGlobalMols_);
158	–	#endif
159	–
162		}
163
164		SimInfo::~SimInfo() {
#	Line 194 \| Line 196 \| namespace oopse {
196		nBonds_ += mol->getNBonds();
197		nBends_ += mol->getNBends();
198		nTorsions_ += mol->getNTorsions();
199	+	nInversions_ += mol->getNInversions();
200		nRigidBodies_ += mol->getNRigidBodies();
201		nIntegrableObjects_ += mol->getNIntegrableObjects();
202		nCutoffGroups_ += mol->getNCutoffGroups();
#	Line 219 \| Line 222 \| namespace oopse {
222		nBonds_ -= mol->getNBonds();
223		nBends_ -= mol->getNBends();
224		nTorsions_ -= mol->getNTorsions();
225	+	nInversions_ -= mol->getNInversions();
226		nRigidBodies_ -= mol->getNRigidBodies();
227		nIntegrableObjects_ -= mol->getNIntegrableObjects();
228		nCutoffGroups_ -= mol->getNCutoffGroups();
#	Line 354 \| Line 358 \| namespace oopse {
358		std::vector<Bond*>::iterator bondIter;
359		std::vector<Bend*>::iterator bendIter;
360		std::vector<Torsion*>::iterator torsionIter;
361	+	std::vector<Inversion*>::iterator inversionIter;
362		Bond* bond;
363		Bend* bend;
364		Torsion* torsion;
365	+	Inversion* inversion;
366		int a;
367		int b;
368		int c;
#	Line 449 \| Line 455 \| namespace oopse {
455		*/
456		}
457
458	+	for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
459	+	inversion = mol->nextInversion(inversionIter)) {
460	+	a = inversion->getAtomA()->getGlobalIndex();
461	+	b = inversion->getAtomB()->getGlobalIndex();
462	+	c = inversion->getAtomC()->getGlobalIndex();
463	+	d = inversion->getAtomD()->getGlobalIndex();
464	+	std::set<int> rigidSetA = getRigidSet(a, atomGroups);
465	+	std::set<int> rigidSetB = getRigidSet(b, atomGroups);
466	+	std::set<int> rigidSetC = getRigidSet(c, atomGroups);
467	+	std::set<int> rigidSetD = getRigidSet(d, atomGroups);
468	+
469	+	exclude_.addPairs(rigidSetA, rigidSetB);
470	+	exclude_.addPairs(rigidSetA, rigidSetC);
471	+	exclude_.addPairs(rigidSetA, rigidSetD);
472	+	exclude_.addPairs(rigidSetB, rigidSetC);
473	+	exclude_.addPairs(rigidSetB, rigidSetD);
474	+	exclude_.addPairs(rigidSetC, rigidSetD);
475	+
476	+	/*
477	+	exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
478	+	exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
479	+	exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
480	+	exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
481	+	exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
482	+	exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
483	+
484	+
485	+	exclude_.addPair(a, b);
486	+	exclude_.addPair(a, c);
487	+	exclude_.addPair(a, d);
488	+	exclude_.addPair(b, c);
489	+	exclude_.addPair(b, d);
490	+	exclude_.addPair(c, d);
491	+	*/
492	+	}
493	+
494		for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
495		std::vector<Atom*> atoms = rb->getAtoms();
496		for (int i = 0; i < atoms.size() -1 ; ++i) {
#	Line 466 \| Line 508 \| namespace oopse {
508		std::vector<Bond*>::iterator bondIter;
509		std::vector<Bend*>::iterator bendIter;
510		std::vector<Torsion*>::iterator torsionIter;
511	+	std::vector<Inversion*>::iterator inversionIter;
512		Bond* bond;
513		Bend* bend;
514		Torsion* torsion;
515	+	Inversion* inversion;
516		int a;
517		int b;
518		int c;
#	Line 562 \| Line 606 \| namespace oopse {
606		*/
607		}
608
609	+	for (inversion= mol->beginInversion(inversionIter); inversion != NULL; inversion = mol->nextInversion(inversionIter)) {
610	+	a = inversion->getAtomA()->getGlobalIndex();
611	+	b = inversion->getAtomB()->getGlobalIndex();
612	+	c = inversion->getAtomC()->getGlobalIndex();
613	+	d = inversion->getAtomD()->getGlobalIndex();
614	+
615	+	std::set<int> rigidSetA = getRigidSet(a, atomGroups);
616	+	std::set<int> rigidSetB = getRigidSet(b, atomGroups);
617	+	std::set<int> rigidSetC = getRigidSet(c, atomGroups);
618	+	std::set<int> rigidSetD = getRigidSet(d, atomGroups);
619	+
620	+	exclude_.removePairs(rigidSetA, rigidSetB);
621	+	exclude_.removePairs(rigidSetA, rigidSetC);
622	+	exclude_.removePairs(rigidSetA, rigidSetD);
623	+	exclude_.removePairs(rigidSetB, rigidSetC);
624	+	exclude_.removePairs(rigidSetB, rigidSetD);
625	+	exclude_.removePairs(rigidSetC, rigidSetD);
626	+
627	+	/*
628	+	exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
629	+	exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
630	+	exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
631	+	exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
632	+	exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
633	+	exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
634	+
635	+
636	+	exclude_.removePair(a, b);
637	+	exclude_.removePair(a, c);
638	+	exclude_.removePair(a, d);
639	+	exclude_.removePair(b, c);
640	+	exclude_.removePair(b, d);
641	+	exclude_.removePair(c, d);
642	+	*/
643	+	}
644	+
645		for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
646		std::vector<Atom*> atoms = rb->getAtoms();
647		for (int i = 0; i < atoms.size() -1 ; ++i) {
#	Line 600 \| Line 680 \| namespace oopse {
680		/** @deprecate */
681		int isError = 0;
682
683	+	setupCutoff();
684	+
685		setupElectrostaticSummationMethod( isError );
686		setupSwitchingFunction();
687	+	setupAccumulateBoxDipole();
688
689		if(isError){
690		sprintf( painCave.errMsg,
#	Line 609 \| Line 692 \| namespace oopse {
692		painCave.isFatal = 1;
693		simError();
694		}
612	–
613	–
614	–	setupCutoff();
695
696		calcNdf();
697		calcNdfRaw();
#	Line 661 \| Line 741 \| namespace oopse {
741		int usePBC = simParams_->getUsePeriodicBoundaryConditions();
742		int useRF;
743		int useSF;
744	+	int useSP;
745	+	int useBoxDipole;
746	+
747		std::string myMethod;
748
749		// set the useRF logical
750		useRF = 0;
751		useSF = 0;
752	+	useSP = 0;
753
754
755		if (simParams_->haveElectrostaticSummationMethod()) {
756		std::string myMethod = simParams_->getElectrostaticSummationMethod();
757		toUpper(myMethod);
758	<	if (myMethod == "REACTION_FIELD") {
759	<	useRF=1;
760	<	} else {
761	<	if (myMethod == "SHIFTED_FORCE") {
762	<	useSF = 1;
763	<	}
758	>	if (myMethod == "REACTION_FIELD"){
759	>	useRF = 1;
760	>	} else if (myMethod == "SHIFTED_FORCE"){
761	>	useSF = 1;
762	>	} else if (myMethod == "SHIFTED_POTENTIAL"){
763	>	useSP = 1;
764		}
765		}
766	+
767	+	if (simParams_->haveAccumulateBoxDipole())
768	+	if (simParams_->getAccumulateBoxDipole())
769	+	useBoxDipole = 1;
770
771	+	useAtomicVirial_ = simParams_->getUseAtomicVirial();
772	+
773		//loop over all of the atom types
774		for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
775		useLennardJones \|= (*i)->isLennardJones();
#	Line 749 \| Line 839 \| namespace oopse {
839		MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
840
841		temp = useSF;
842	<	MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
842	>	MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
843
844	+	temp = useSP;
845	+	MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
846	+
847	+	temp = useBoxDipole;
848	+	MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
849	+
850	+	temp = useAtomicVirial_;
851	+	MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
852	+
853		#endif
854
855		fInfo_.SIM_uses_PBC = usePBC;
#	Line 768 \| Line 867 \| namespace oopse {
867		fInfo_.SIM_uses_FLARB = useFLARB;
868		fInfo_.SIM_uses_RF = useRF;
869		fInfo_.SIM_uses_SF = useSF;
870	<
871	<	if( myMethod == "REACTION_FIELD") {
872	<
774	<	if (simParams_->haveDielectric()) {
775	<	fInfo_.dielect = simParams_->getDielectric();
776	<	} else {
777	<	sprintf(painCave.errMsg,
778	<	"SimSetup Error: No Dielectric constant was set.\n"
779	<	"\tYou are trying to use Reaction Field without"
780	<	"\tsetting a dielectric constant!\n");
781	<	painCave.isFatal = 1;
782	<	simError();
783	<	}
784	<	}
785	<
870	>	fInfo_.SIM_uses_SP = useSP;
871	>	fInfo_.SIM_uses_BoxDipole = useBoxDipole;
872	>	fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_;
873		}
874
875		void SimInfo::setupFortranSim() {
#	Line 799 \| Line 886 \| namespace oopse {
886		}
887
888		//calculate mass ratio of cutoff group
889	<	std::vector<double> mfact;
889	>	std::vector<RealType> mfact;
890		SimInfo::MoleculeIterator mi;
891		Molecule* mol;
892		Molecule::CutoffGroupIterator ci;
893		CutoffGroup* cg;
894		Molecule::AtomIterator ai;
895		Atom* atom;
896	<	double totalMass;
896	>	RealType totalMass;
897
898		//to avoid memory reallocation, reserve enough space for mfact
899		mfact.reserve(getNCutoffGroups());
#	Line 849 \| Line 936 \| namespace oopse {
936		int nGlobalExcludes = 0;
937		int* globalExcludes = NULL;
938		int* excludeList = exclude_.getExcludeList();
939	<	setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
940	<	&nGlobalExcludes, globalExcludes, &molMembershipArray[0],
941	<	&mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
942	<
939	>	setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0],
940	>	&nExclude, excludeList , &nGlobalExcludes, globalExcludes,
941	>	&molMembershipArray[0], &mfact[0], &nCutoffGroups_,
942	>	&fortranGlobalGroupMembership[0], &isError);
943	>
944		if( isError ){
945	<
945	>
946		sprintf( painCave.errMsg,
947		"There was an error setting the simulation information in fortran.\n" );
948		painCave.isFatal = 1;
949		painCave.severity = OOPSE_ERROR;
950		simError();
951		}
952	<
953	<	#ifdef IS_MPI
952	>
953	>
954		sprintf( checkPointMsg,
955		"succesfully sent the simulation information to fortran.\n");
956	<	MPIcheckPoint();
957	<	#endif // is_mpi
956	>
957	>	errorCheckPoint();
958	>
959	>	// Setup number of neighbors in neighbor list if present
960	>	if (simParams_->haveNeighborListNeighbors()) {
961	>	int nlistNeighbors = simParams_->getNeighborListNeighbors();
962	>	setNeighbors(&nlistNeighbors);
963	>	}
964	>
965	>
966		}
967
968
873	–	#ifdef IS_MPI
969		void SimInfo::setupFortranParallel() {
970	<
970	>	#ifdef IS_MPI
971		//SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
972		std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
973		std::vector<int> localToGlobalCutoffGroupIndex;
#	Line 922 \| Line 1017 \| namespace oopse {
1017		}
1018
1019		sprintf(checkPointMsg, " mpiRefresh successful.\n");
1020	<	MPIcheckPoint();
1020	>	errorCheckPoint();
1021
1022	<
1022	>	#endif
1023		}
1024
930	–	#endif
931	–
1025		void SimInfo::setupCutoff() {
1026
1027		ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
#	Line 936 \| Line 1029 \| namespace oopse {
1029		// Check the cutoff policy
1030		int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default
1031
1032	+	// Set LJ shifting bools to false
1033	+	ljsp_ = false;
1034	+	ljsf_ = false;
1035	+
1036		std::string myPolicy;
1037		if (forceFieldOptions_.haveCutoffPolicy()){
1038		myPolicy = forceFieldOptions_.getCutoffPolicy();
#	Line 966 \| Line 1063 \| namespace oopse {
1063		notifyFortranCutoffPolicy(&cp);
1064
1065		// Check the Skin Thickness for neighborlists
1066	<	double skin;
1066	>	RealType skin;
1067		if (simParams_->haveSkinThickness()) {
1068		skin = simParams_->getSkinThickness();
1069		notifyFortranSkinThickness(&skin);
#	Line 999 \| Line 1096 \| namespace oopse {
1096		simError();
1097		}
1098		}
1099	+
1100	+	if (simParams_->haveElectrostaticSummationMethod()) {
1101	+	std::string myMethod = simParams_->getElectrostaticSummationMethod();
1102	+	toUpper(myMethod);
1103	+
1104	+	if (myMethod == "SHIFTED_POTENTIAL") {
1105	+	ljsp_ = true;
1106	+	} else if (myMethod == "SHIFTED_FORCE") {
1107	+	ljsf_ = true;
1108	+	}
1109	+	}
1110	+	notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
1111
1003	–	notifyFortranCutoffs(&rcut_, &rsw_);
1004	–
1112		} else {
1113
1114		// For electrostatic atoms, we'll assume a large safe value:
#	Line 1017 \| Line 1124 \| namespace oopse {
1124		if (simParams_->haveElectrostaticSummationMethod()) {
1125		std::string myMethod = simParams_->getElectrostaticSummationMethod();
1126		toUpper(myMethod);
1127	<	if (myMethod == "SHIFTED_POTENTIAL" \|\| myMethod == "SHIFTED_FORCE") {
1127	>
1128	>	// For the time being, we're tethering the LJ shifted behavior to the
1129	>	// electrostaticSummationMethod keyword options
1130	>	if (myMethod == "SHIFTED_POTENTIAL") {
1131	>	ljsp_ = true;
1132	>	} else if (myMethod == "SHIFTED_FORCE") {
1133	>	ljsf_ = true;
1134	>	}
1135	>	if (myMethod == "SHIFTED_POTENTIAL" \|\| myMethod == "SHIFTED_FORCE") {
1136		if (simParams_->haveSwitchingRadius()){
1137		sprintf(painCave.errMsg,
1138		"SimInfo Warning: A value was set for the switchingRadius\n"
#	Line 1040 \| Line 1155 \| namespace oopse {
1155		simError();
1156		rsw_ = 0.85 * rcut_;
1157		}
1158	<	notifyFortranCutoffs(&rcut_, &rsw_);
1158	>
1159	>	notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
1160	>
1161		} else {
1162		// We didn't set rcut explicitly, and we don't have electrostatic atoms, so
1163		// We'll punt and let fortran figure out the cutoffs later.
#	Line 1056 \| Line 1173 \| namespace oopse {
1173		int errorOut;
1174		int esm = NONE;
1175		int sm = UNDAMPED;
1176	<	double alphaVal;
1177	<	double dielectric;
1178	<
1176	>	RealType alphaVal;
1177	>	RealType dielectric;
1178	>
1179		errorOut = isError;
1063	–	alphaVal = simParams_->getDampingAlpha();
1064	–	dielectric = simParams_->getDielectric();
1180
1181		if (simParams_->haveElectrostaticSummationMethod()) {
1182		std::string myMethod = simParams_->getElectrostaticSummationMethod();
#	Line 1078 \| Line 1193 \| namespace oopse {
1193		if (myMethod == "SHIFTED_FORCE") {
1194		esm = SHIFTED_FORCE;
1195		} else {
1196	<	if (myMethod == "REACTION_FIELD") {
1196	>	if (myMethod == "REACTION_FIELD") {
1197		esm = REACTION_FIELD;
1198	+	dielectric = simParams_->getDielectric();
1199	+	if (!simParams_->haveDielectric()) {
1200	+	// throw warning
1201	+	sprintf( painCave.errMsg,
1202	+	"SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
1203	+	"\tA default value of %f will be used for the dielectric.\n", dielectric);
1204	+	painCave.isFatal = 0;
1205	+	simError();
1206	+	}
1207		} else {
1208		// throw error
1209		sprintf( painCave.errMsg,
#	Line 1106 \| Line 1230 \| namespace oopse {
1230		if (myScreen == "DAMPED") {
1231		sm = DAMPED;
1232		if (!simParams_->haveDampingAlpha()) {
1233	<	//throw error
1233	>	// first set a cutoff dependent alpha value
1234	>	// we assume alpha depends linearly with rcut from 0 to 20.5 ang
1235	>	alphaVal = 0.5125 - rcut_* 0.025;
1236	>	// for values rcut > 20.5, alpha is zero
1237	>	if (alphaVal < 0) alphaVal = 0;
1238	>
1239	>	// throw warning
1240		sprintf( painCave.errMsg,
1241		"SimInfo warning: dampingAlpha was not specified in the input file.\n"
1242	<	"\tA default value of %f (1/ang) will be used.\n", alphaVal);
1242	>	"\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_);
1243		painCave.isFatal = 0;
1244		simError();
1245	+	} else {
1246	+	alphaVal = simParams_->getDampingAlpha();
1247		}
1248	+
1249		} else {
1250		// throw error
1251		sprintf( painCave.errMsg,
#	Line 1161 \| Line 1294 \| namespace oopse {
1294
1295		}
1296
1297	+	void SimInfo::setupAccumulateBoxDipole() {
1298	+
1299	+	// we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
1300	+	if ( simParams_->haveAccumulateBoxDipole() )
1301	+	if ( simParams_->getAccumulateBoxDipole() ) {
1302	+	setAccumulateBoxDipole();
1303	+	calcBoxDipole_ = true;
1304	+	}
1305	+
1306	+	}
1307	+
1308		void SimInfo::addProperty(GenericData* genData) {
1309		properties_.addProperty(genData);
1310		}
#	Line 1217 \| Line 1361 \| namespace oopse {
1361		Molecule* mol;
1362
1363		Vector3d comVel(0.0);
1364	<	double totalMass = 0.0;
1364	>	RealType totalMass = 0.0;
1365
1366
1367		for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1368	<	double mass = mol->getMass();
1368	>	RealType mass = mol->getMass();
1369		totalMass += mass;
1370		comVel += mass * mol->getComVel();
1371		}
1372
1373		#ifdef IS_MPI
1374	<	double tmpMass = totalMass;
1374	>	RealType tmpMass = totalMass;
1375		Vector3d tmpComVel(comVel);
1376	<	MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1377	<	MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1376	>	MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1377	>	MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1378		#endif
1379
1380		comVel /= totalMass;
#	Line 1243 \| Line 1387 \| namespace oopse {
1387		Molecule* mol;
1388
1389		Vector3d com(0.0);
1390	<	double totalMass = 0.0;
1390	>	RealType totalMass = 0.0;
1391
1392		for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1393	<	double mass = mol->getMass();
1393	>	RealType mass = mol->getMass();
1394		totalMass += mass;
1395		com += mass * mol->getCom();
1396		}
1397
1398		#ifdef IS_MPI
1399	<	double tmpMass = totalMass;
1399	>	RealType tmpMass = totalMass;
1400		Vector3d tmpCom(com);
1401	<	MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1402	<	MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1401	>	MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1402	>	MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1403		#endif
1404
1405		com /= totalMass;
#	Line 1279 \| Line 1423 \| namespace oopse {
1423		Molecule* mol;
1424
1425
1426	<	double totalMass = 0.0;
1426	>	RealType totalMass = 0.0;
1427
1428
1429		for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1430	<	double mass = mol->getMass();
1430	>	RealType mass = mol->getMass();
1431		totalMass += mass;
1432		com += mass * mol->getCom();
1433		comVel += mass * mol->getComVel();
1434		}
1435
1436		#ifdef IS_MPI
1437	<	double tmpMass = totalMass;
1437	>	RealType tmpMass = totalMass;
1438		Vector3d tmpCom(com);
1439		Vector3d tmpComVel(comVel);
1440	<	MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1441	<	MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1442	<	MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1440	>	MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1441	>	MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1442	>	MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1443		#endif
1444
1445		com /= totalMass;
#	Line 1314 \| Line 1458 \| namespace oopse {
1458		void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1459
1460
1461	<	double xx = 0.0;
1462	<	double yy = 0.0;
1463	<	double zz = 0.0;
1464	<	double xy = 0.0;
1465	<	double xz = 0.0;
1466	<	double yz = 0.0;
1461	>	RealType xx = 0.0;
1462	>	RealType yy = 0.0;
1463	>	RealType zz = 0.0;
1464	>	RealType xy = 0.0;
1465	>	RealType xz = 0.0;
1466	>	RealType yz = 0.0;
1467		Vector3d com(0.0);
1468		Vector3d comVel(0.0);
1469
#	Line 1331 \| Line 1475 \| namespace oopse {
1475		Vector3d thisq(0.0);
1476		Vector3d thisv(0.0);
1477
1478	<	double thisMass = 0.0;
1478	>	RealType thisMass = 0.0;
1479
1480
1481
#	Line 1369 \| Line 1513 \| namespace oopse {
1513		#ifdef IS_MPI
1514		Mat3x3d tmpI(inertiaTensor);
1515		Vector3d tmpAngMom;
1516	<	MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1517	<	MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1516	>	MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1517	>	MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1518		#endif
1519
1520		return;
#	Line 1391 \| Line 1535 \| namespace oopse {
1535		Vector3d thisr(0.0);
1536		Vector3d thisp(0.0);
1537
1538	<	double thisMass;
1538	>	RealType thisMass;
1539
1540		for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1541		thisMass = mol->getMass();
#	Line 1404 \| Line 1548 \| namespace oopse {
1548
1549		#ifdef IS_MPI
1550		Vector3d tmpAngMom;
1551	<	MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1551	>	MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1552		#endif
1553
1554		return angularMomentum;
1555		}
1556
1557	<
1557	>	StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1558	>	return IOIndexToIntegrableObject.at(index);
1559	>	}
1560	>
1561	>	void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
1562	>	IOIndexToIntegrableObject= v;
1563	>	}
1564	>
1565	>	/* Returns the Volume of the simulation based on a ellipsoid with semi-axes
1566	>	based on the radius of gyration V=4/3PiR_1R_2R_3
1567	>	where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to
1568	>	V = 4/3Pi(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
1569	>	*/
1570	>	void SimInfo::getGyrationalVolume(RealType &volume){
1571	>	Mat3x3d intTensor;
1572	>	RealType det;
1573	>	Vector3d dummyAngMom;
1574	>	RealType sysconstants;
1575	>	RealType geomCnst;
1576	>
1577	>	geomCnst = 3.0/2.0;
1578	>	/* Get the inertial tensor and angular momentum for free*/
1579	>	getInertiaTensor(intTensor,dummyAngMom);
1580	>
1581	>	det = intTensor.determinant();
1582	>	sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1583	>	volume = 4.0/3.0NumericConstant::PIpow(sysconstants,3.0/2.0)*sqrt(det);
1584	>	return;
1585	>	}
1586	>
1587	>	void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
1588	>	Mat3x3d intTensor;
1589	>	Vector3d dummyAngMom;
1590	>	RealType sysconstants;
1591	>	RealType geomCnst;
1592	>
1593	>	geomCnst = 3.0/2.0;
1594	>	/* Get the inertial tensor and angular momentum for free*/
1595	>	getInertiaTensor(intTensor,dummyAngMom);
1596	>
1597	>	detI = intTensor.determinant();
1598	>	sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1599	>	volume = 4.0/3.0NumericConstant::PIpow(sysconstants,3.0/2.0)*sqrt(detI);
1600	>	return;
1601	>	}
1602	>	/*
1603	>	void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
1604	>	assert( v.size() == nAtoms_ + nRigidBodies_);
1605	>	sdByGlobalIndex_ = v;
1606	>	}
1607	>
1608	>	StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
1609	>	//assert(index < nAtoms_ + nRigidBodies_);
1610	>	return sdByGlobalIndex_.at(index);
1611	>	}
1612	>	*/
1613		}//end namespace oopse
1614

Diff Legend

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

Comparing trunk/src/brains/SimInfo.cpp (file contents): Revision 945 by gezelter, Tue Apr 25 02:09:01 2006 UTC vs. Revision 1277 by gezelter, Mon Jul 14 12:35:58 2008 UTC

Diff Legend

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 945 by gezelter, Tue Apr 25 02:09:01 2006 UTC vs.
Revision 1277 by gezelter, Mon Jul 14 12:35:58 2008 UTC