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

Comparing branches/development/src/brains/SimInfo.cpp (file contents):
Revision 1532 by gezelter, Wed Dec 29 19:59:21 2010 UTC vs.
Revision 1725 by gezelter, Sat May 26 18:13:43 2012 UTC

#	Line 36 \| Line 36
36		* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37		* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38		* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	<	* [4] Vardeman & Gezelter, in progress (2009).
39	>	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	>	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43		/**
#	Line 54 \| Line 55
55		#include "math/Vector3.hpp"
56		#include "primitives/Molecule.hpp"
57		#include "primitives/StuntDouble.hpp"
57	–	#include "UseTheForce/doForces_interface.h"
58	–	#include "UseTheForce/DarkSide/neighborLists_interface.h"
58		#include "utils/MemoryUtils.hpp"
59		#include "utils/simError.h"
60		#include "selection/SelectionManager.hpp"
61		#include "io/ForceFieldOptions.hpp"
62	<	#include "UseTheForce/ForceField.hpp"
62	>	#include "brains/ForceField.hpp"
63		#include "nonbonded/SwitchingFunction.hpp"
65	–
66	–
64		#ifdef IS_MPI
65	<	#include "UseTheForce/mpiComponentPlan.h"
66	<	#include "UseTheForce/DarkSide/simParallel_interface.h"
70	<	#endif
65	>	#include <mpi.h>
66	>	#endif
67
68		using namespace std;
69		namespace OpenMD {
#	Line 76 \| Line 72 \| namespace OpenMD {
72		forceField_(ff), simParams_(simParams),
73		ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
74		nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
75	<	nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
75	>	nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0),
76		nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0),
77		nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0),
78	<	nConstraints_(0), sman_(NULL), fortranInitialized_(false),
78	>	nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false),
79		calcBoxDipole_(false), useAtomicVirial_(true) {
80
81		MoleculeStamp* molStamp;
#	Line 133 \| Line 129 \| namespace OpenMD {
129		//equal to the total number of atoms minus number of atoms belong to
130		//cutoff group defined in meta-data file plus the number of cutoff
131		//groups defined in meta-data file
132	+
133		nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
134
135		//every free atom (atom does not belong to rigid bodies) is an
#	Line 228 \| Line 225 \| namespace OpenMD {
225
226
227		void SimInfo::calcNdf() {
228	<	int ndf_local;
228	>	int ndf_local, nfq_local;
229		MoleculeIterator i;
230		vector<StuntDouble*>::iterator j;
231	+	vector<Atom*>::iterator k;
232	+
233		Molecule* mol;
234		StuntDouble* integrableObject;
235	+	Atom* atom;
236
237		ndf_local = 0;
238	+	nfq_local = 0;
239
240		for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
241		for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
#	Line 249 \| Line 250 \| namespace OpenMD {
250		ndf_local += 3;
251		}
252		}
252	–
253		}
254	+	for (atom = mol->beginFluctuatingCharge(k); atom != NULL;
255	+	atom = mol->nextFluctuatingCharge(k)) {
256	+	if (atom->isFluctuatingCharge()) {
257	+	nfq_local++;
258	+	}
259	+	}
260		}
261
262		// n_constraints is local, so subtract them on each processor
#	Line 258 \| Line 264 \| namespace OpenMD {
264
265		#ifdef IS_MPI
266		MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
267	+	MPI_Allreduce(&nfq_local,&nGlobalFluctuatingCharges_,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
268		#else
269		ndf_ = ndf_local;
270	+	nGlobalFluctuatingCharges_ = nfq_local;
271		#endif
272
273		// nZconstraints_ is global, as are the 3 COM translations for the
#	Line 275 \| Line 283 \| namespace OpenMD {
283		fdf_ = fdf_local;
284		#endif
285		return fdf_;
286	+	}
287	+
288	+	unsigned int SimInfo::getNLocalCutoffGroups(){
289	+	int nLocalCutoffAtoms = 0;
290	+	Molecule* mol;
291	+	MoleculeIterator mi;
292	+	CutoffGroup* cg;
293	+	Molecule::CutoffGroupIterator ci;
294	+
295	+	for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
296	+
297	+	for (cg = mol->beginCutoffGroup(ci); cg != NULL;
298	+	cg = mol->nextCutoffGroup(ci)) {
299	+	nLocalCutoffAtoms += cg->getNumAtom();
300	+
301	+	}
302	+	}
303	+
304	+	return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_;
305		}
306
307		void SimInfo::calcNdfRaw() {
#	Line 657 \| Line 684 \| namespace OpenMD {
684		/**
685		* update
686		*
687	<	* Performs the global checks and variable settings after the objects have been
688	<	* created.
687	>	* Performs the global checks and variable settings after the
688	>	* objects have been created.
689		*
690		*/
691	<	void SimInfo::update() {
665	<
691	>	void SimInfo::update() {
692		setupSimVariables();
667	–	setupCutoffs();
668	–	setupSwitching();
669	–	setupElectrostatics();
670	–	setupNeighborlists();
671	–
672	–	#ifdef IS_MPI
673	–	setupFortranParallel();
674	–	#endif
675	–	setupFortranSim();
676	–	fortranInitialized_ = true;
677	–
693		calcNdf();
694		calcNdfRaw();
695		calcNdfTrans();
696		}
697
698	+	/**
699	+	* getSimulatedAtomTypes
700	+	*
701	+	* Returns an STL set of AtomType* that are actually present in this
702	+	* simulation. Must query all processors to assemble this information.
703	+	*
704	+	*/
705		set<AtomType*> SimInfo::getSimulatedAtomTypes() {
706		SimInfo::MoleculeIterator mi;
707		Molecule* mol;
#	Line 687 \| Line 709 \| namespace OpenMD {
709		Atom* atom;
710		set<AtomType*> atomTypes;
711
712	<	for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
713	<	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
712	>	for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
713	>	for(atom = mol->beginAtom(ai); atom != NULL;
714	>	atom = mol->nextAtom(ai)) {
715		atomTypes.insert(atom->getAtomType());
716		}
717		}
695	–	return atomTypes;
696	–	}
697	–
698	–	/**
699	–	* setupCutoffs
700	–	*
701	–	* Sets the values of cutoffRadius and cutoffMethod
702	–	*
703	–	* cutoffRadius : realType
704	–	* If the cutoffRadius was explicitly set, use that value.
705	–	* If the cutoffRadius was not explicitly set:
706	–	* Are there electrostatic atoms? Use 12.0 Angstroms.
707	–	* No electrostatic atoms? Poll the atom types present in the
708	–	* simulation for suggested cutoff values (e.g. 2.5 * sigma).
709	–	* Use the maximum suggested value that was found.
710	–	*
711	–	* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
712	–	* If cutoffMethod was explicitly set, use that choice.
713	–	* If cutoffMethod was not explicitly set, use SHIFTED_FORCE
714	–	*/
715	–	void SimInfo::setupCutoffs() {
718
719	<	if (simParams_->haveCutoffRadius()) {
718	<	cutoffRadius_ = simParams_->getCutoffRadius();
719	<	} else {
720	<	if (usesElectrostaticAtoms_) {
721	<	sprintf(painCave.errMsg,
722	<	"SimInfo: No value was set for the cutoffRadius.\n"
723	<	"\tOpenMD will use a default value of 12.0 angstroms"
724	<	"\tfor the cutoffRadius.\n");
725	<	painCave.isFatal = 0;
726	<	painCave.severity = OPENMD_INFO;
727	<	simError();
728	<	cutoffRadius_ = 12.0;
729	<	} else {
730	<	RealType thisCut;
731	<	set<AtomType*>::iterator i;
732	<	set<AtomType*> atomTypes;
733	<	atomTypes = getSimulatedAtomTypes();
734	<	for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
735	<	thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i));
736	<	cutoffRadius_ = max(thisCut, cutoffRadius_);
737	<	}
738	<	sprintf(painCave.errMsg,
739	<	"SimInfo: No value was set for the cutoffRadius.\n"
740	<	"\tOpenMD will use %lf angstroms.\n",
741	<	cutoffRadius_);
742	<	painCave.isFatal = 0;
743	<	painCave.severity = OPENMD_INFO;
744	<	simError();
745	<	}
746	<	}
719	>	#ifdef IS_MPI
720
721	<	InteractionManager::Instance()->setCutoffRadius(cutoffRadius_);
722	<
723	<	map<string, CutoffMethod> stringToCutoffMethod;
724	<	stringToCutoffMethod["HARD"] = HARD;
725	<	stringToCutoffMethod["SWITCHING_FUNCTION"] = SWITCHING_FUNCTION;
726	<	stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;
727	<	stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
721	>	// loop over the found atom types on this processor, and add their
722	>	// numerical idents to a vector:
723	>
724	>	vector<int> foundTypes;
725	>	set<AtomType*>::iterator i;
726	>	for (i = atomTypes.begin(); i != atomTypes.end(); ++i)
727	>	foundTypes.push_back( (*i)->getIdent() );
728	>
729	>	// count_local holds the number of found types on this processor
730	>	int count_local = foundTypes.size();
731	>
732	>	int nproc = MPI::COMM_WORLD.Get_size();
733	>
734	>	// we need arrays to hold the counts and displacement vectors for
735	>	// all processors
736	>	vector<int> counts(nproc, 0);
737	>	vector<int> disps(nproc, 0);
738	>
739	>	// fill the counts array
740	>	MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0],
741	>	1, MPI::INT);
742
743	<	if (simParams_->haveCutoffMethod()) {
744	<	string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
745	<	map<string, CutoffMethod>::iterator i;
746	<	i = stringToCutoffMethod.find(cutMeth);
747	<	if (i == stringToCutoffMethod.end()) {
748	<	sprintf(painCave.errMsg,
762	<	"SimInfo: Could not find chosen cutoffMethod %s\n"
763	<	"\tShould be one of: "
764	<	"HARD, SWITCHING_FUNCTION, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
765	<	cutMeth.c_str());
766	<	painCave.isFatal = 1;
767	<	painCave.severity = OPENMD_ERROR;
768	<	simError();
769	<	} else {
770	<	cutoffMethod_ = i->second;
771	<	}
772	<	} else {
773	<	sprintf(painCave.errMsg,
774	<	"SimInfo: No value was set for the cutoffMethod.\n"
775	<	"\tOpenMD will use SHIFTED_FORCE.\n");
776	<	painCave.isFatal = 0;
777	<	painCave.severity = OPENMD_INFO;
778	<	simError();
779	<	cutoffMethod_ = SHIFTED_FORCE;
743	>	// use the processor counts to compute the displacement array
744	>	disps[0] = 0;
745	>	int totalCount = counts[0];
746	>	for (int iproc = 1; iproc < nproc; iproc++) {
747	>	disps[iproc] = disps[iproc-1] + counts[iproc-1];
748	>	totalCount += counts[iproc];
749		}
750
751	<	InteractionManager::Instance()->setCutoffMethod(cutoffMethod_);
752	<	}
784	<
785	<	/**
786	<	* setupSwitching
787	<	*
788	<	* Sets the values of switchingRadius and
789	<	* If the switchingRadius was explicitly set, use that value (but check it)
790	<	* If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
791	<	*/
792	<	void SimInfo::setupSwitching() {
751	>	// we need a (possibly redundant) set of all found types:
752	>	vector<int> ftGlobal(totalCount);
753
754	<	if (simParams_->haveSwitchingRadius()) {
755	<	switchingRadius_ = simParams_->getSwitchingRadius();
756	<	if (switchingRadius_ > cutoffRadius_) {
757	<	sprintf(painCave.errMsg,
798	<	"SimInfo: switchingRadius (%f) is larger than cutoffRadius(%f)\n",
799	<	switchingRadius_, cutoffRadius_);
800	<	painCave.isFatal = 1;
801	<	painCave.severity = OPENMD_ERROR;
802	<	simError();
803	<	}
804	<	} else {
805	<	switchingRadius_ = 0.85 * cutoffRadius_;
806	<	sprintf(painCave.errMsg,
807	<	"SimInfo: No value was set for the switchingRadius.\n"
808	<	"\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
809	<	"\tswitchingRadius = %f. for this simulation\n", switchingRadius_);
810	<	painCave.isFatal = 0;
811	<	painCave.severity = OPENMD_WARNING;
812	<	simError();
813	<	}
814	<
815	<	InteractionManager::Instance()->setSwitchingRadius(switchingRadius_);
754	>	// now spray out the foundTypes to all the other processors:
755	>	MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT,
756	>	&ftGlobal[0], &counts[0], &disps[0],
757	>	MPI::INT);
758
759	<	SwitchingFunctionType ft;
818	<
819	<	if (simParams_->haveSwitchingFunctionType()) {
820	<	string funcType = simParams_->getSwitchingFunctionType();
821	<	toUpper(funcType);
822	<	if (funcType == "CUBIC") {
823	<	ft = cubic;
824	<	} else {
825	<	if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
826	<	ft = fifth_order_poly;
827	<	} else {
828	<	// throw error
829	<	sprintf( painCave.errMsg,
830	<	"SimInfo : Unknown switchingFunctionType. (Input file specified %s .)\n"
831	<	"\tswitchingFunctionType must be one of: "
832	<	"\"cubic\" or \"fifth_order_polynomial\".",
833	<	funcType.c_str() );
834	<	painCave.isFatal = 1;
835	<	painCave.severity = OPENMD_ERROR;
836	<	simError();
837	<	}
838	<	}
839	<	}
759	>	vector<int>::iterator j;
760
761	<	InteractionManager::Instance()->setSwitchingFunctionType(ft);
762	<	}
761	>	// foundIdents is a stl set, so inserting an already found ident
762	>	// will have no effect.
763	>	set<int> foundIdents;
764
765	<	/**
766	<	* setupSkinThickness
767	<	*
768	<	* If the skinThickness was explicitly set, use that value (but check it)
769	<	* If the skinThickness was not explicitly set: use 1.0 angstroms
770	<	*/
771	<	void SimInfo::setupSkinThickness() {
772	<	if (simParams_->haveSkinThickness()) {
773	<	skinThickness_ = simParams_->getSkinThickness();
774	<	} else {
775	<	skinThickness_ = 1.0;
776	<	sprintf(painCave.errMsg,
856	<	"SimInfo Warning: No value was set for the skinThickness.\n"
857	<	"\tOpenMD will use a default value of %f Angstroms\n"
858	<	"\tfor this simulation\n", skinThickness_);
859	<	painCave.isFatal = 0;
860	<	simError();
861	<	}
765	>	for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j)
766	>	foundIdents.insert((*j));
767	>
768	>	// now iterate over the foundIdents and get the actual atom types
769	>	// that correspond to these:
770	>	set<int>::iterator it;
771	>	for (it = foundIdents.begin(); it != foundIdents.end(); ++it)
772	>	atomTypes.insert( forceField_->getAtomType((*it)) );
773	>
774	>	#endif
775	>
776	>	return atomTypes;
777		}
778
779	<	void SimInfo::setupSimType() {
779	>	void SimInfo::setupSimVariables() {
780	>	useAtomicVirial_ = simParams_->getUseAtomicVirial();
781	>	// we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
782	>	calcBoxDipole_ = false;
783	>	if ( simParams_->haveAccumulateBoxDipole() )
784	>	if ( simParams_->getAccumulateBoxDipole() ) {
785	>	calcBoxDipole_ = true;
786	>	}
787	>
788		set<AtomType*>::iterator i;
789		set<AtomType*> atomTypes;
790	<	atomTypes = getSimulatedAtomTypes();
868	<
869	<	useAtomicVirial_ = simParams_->getUseAtomicVirial();
870	<
790	>	atomTypes = getSimulatedAtomTypes();
791		int usesElectrostatic = 0;
792		int usesMetallic = 0;
793		int usesDirectional = 0;
794	+	int usesFluctuatingCharges = 0;
795		//loop over all of the atom types
796		for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
797		usesElectrostatic \|= (*i)->isElectrostatic();
798		usesMetallic \|= (*i)->isMetal();
799		usesDirectional \|= (*i)->isDirectional();
800	+	usesFluctuatingCharges \|= (*i)->isFluctuatingCharge();
801		}
802	<
802	>
803		#ifdef IS_MPI
804		int temp;
805		temp = usesDirectional;
806		MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
807	<
807	>
808		temp = usesMetallic;
809		MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
810	<
810	>
811		temp = usesElectrostatic;
812		MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
813	+
814	+	temp = usesFluctuatingCharges;
815	+	MPI_Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
816	+	#else
817	+
818	+	usesDirectionalAtoms_ = usesDirectional;
819	+	usesMetallicAtoms_ = usesMetallic;
820	+	usesElectrostaticAtoms_ = usesElectrostatic;
821	+	usesFluctuatingCharges_ = usesFluctuatingCharges;
822	+
823		#endif
824	<	fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_;
825	<	fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_;
826	<	fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_;
827	<	fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_;
896	<	fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_;
897	<	fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_;
824	>
825	>	requiresPrepair_ = usesMetallicAtoms_ ? true : false;
826	>	requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false;
827	>	requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false;
828		}
829
830	<	void SimInfo::setupFortranSim() {
831	<	int isError;
832	<	int nExclude, nOneTwo, nOneThree, nOneFour;
833	<	vector<int> fortranGlobalGroupMembership;
830	>
831	>	vector<int> SimInfo::getGlobalAtomIndices() {
832	>	SimInfo::MoleculeIterator mi;
833	>	Molecule* mol;
834	>	Molecule::AtomIterator ai;
835	>	Atom* atom;
836	>
837	>	vector<int> GlobalAtomIndices(getNAtoms(), 0);
838
839	<	notifyFortranSkinThickness(&skinThickness_);
839	>	for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
840	>
841	>	for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
842	>	GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex();
843	>	}
844	>	}
845	>	return GlobalAtomIndices;
846	>	}
847
907	–	int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0;
908	–	int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0;
909	–	notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf);
848
849	<	isError = 0;
849	>	vector<int> SimInfo::getGlobalGroupIndices() {
850	>	SimInfo::MoleculeIterator mi;
851	>	Molecule* mol;
852	>	Molecule::CutoffGroupIterator ci;
853	>	CutoffGroup* cg;
854
855	<	//globalGroupMembership_ is filled by SimCreator
856	<	for (int i = 0; i < nGlobalAtoms_; i++) {
857	<	fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
855	>	vector<int> GlobalGroupIndices;
856	>
857	>	for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
858	>
859	>	//local index of cutoff group is trivial, it only depends on the
860	>	//order of travesing
861	>	for (cg = mol->beginCutoffGroup(ci); cg != NULL;
862	>	cg = mol->nextCutoffGroup(ci)) {
863	>	GlobalGroupIndices.push_back(cg->getGlobalIndex());
864	>	}
865		}
866	+	return GlobalGroupIndices;
867	+	}
868
869	+
870	+	void SimInfo::prepareTopology() {
871	+	int nExclude, nOneTwo, nOneThree, nOneFour;
872	+
873		//calculate mass ratio of cutoff group
919	–	vector<RealType> mfact;
874		SimInfo::MoleculeIterator mi;
875		Molecule* mol;
876		Molecule::CutoffGroupIterator ci;
#	Line 925 \| Line 879 \| namespace OpenMD {
879		Atom* atom;
880		RealType totalMass;
881
882	<	//to avoid memory reallocation, reserve enough space for mfact
883	<	mfact.reserve(getNCutoffGroups());
882	>	/**
883	>	* The mass factor is the relative mass of an atom to the total
884	>	* mass of the cutoff group it belongs to. By default, all atoms
885	>	* are their own cutoff groups, and therefore have mass factors of
886	>	* 1. We need some special handling for massless atoms, which
887	>	* will be treated as carrying the entire mass of the cutoff
888	>	* group.
889	>	*/
890	>	massFactors_.clear();
891	>	massFactors_.resize(getNAtoms(), 1.0);
892
893		for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
894	<	for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
894	>	for (cg = mol->beginCutoffGroup(ci); cg != NULL;
895	>	cg = mol->nextCutoffGroup(ci)) {
896
897		totalMass = cg->getMass();
898		for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
899		// Check for massless groups - set mfact to 1 if true
900	<	if (totalMass != 0)
901	<	mfact.push_back(atom->getMass()/totalMass);
900	>	if (totalMass != 0)
901	>	massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass;
902		else
903	<	mfact.push_back( 1.0 );
903	>	massFactors_[atom->getLocalIndex()] = 1.0;
904		}
905		}
906		}
907
908	<	//fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
946	<	vector<int> identArray;
908	>	// Build the identArray_
909
910	<	//to avoid memory reallocation, reserve enough space identArray
911	<	identArray.reserve(getNAtoms());
950	<
910	>	identArray_.clear();
911	>	identArray_.reserve(getNAtoms());
912		for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
913		for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
914	<	identArray.push_back(atom->getIdent());
914	>	identArray_.push_back(atom->getIdent());
915		}
916		}
956	–
957	–	//fill molMembershipArray
958	–	//molMembershipArray is filled by SimCreator
959	–	vector<int> molMembershipArray(nGlobalAtoms_);
960	–	for (int i = 0; i < nGlobalAtoms_; i++) {
961	–	molMembershipArray[i] = globalMolMembership_[i] + 1;
962	–	}
917
918	<	//setup fortran simulation
918	>	//scan topology
919
920		nExclude = excludedInteractions_.getSize();
921		nOneTwo = oneTwoInteractions_.getSize();
#	Line 973 \| Line 927 \| namespace OpenMD {
927		int* oneThreeList = oneThreeInteractions_.getPairList();
928		int* oneFourList = oneFourInteractions_.getPairList();
929
930	<	setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0],
977	<	&nExclude, excludeList,
978	<	&nOneTwo, oneTwoList,
979	<	&nOneThree, oneThreeList,
980	<	&nOneFour, oneFourList,
981	<	&molMembershipArray[0], &mfact[0], &nCutoffGroups_,
982	<	&fortranGlobalGroupMembership[0], &isError);
983	<
984	<	if( isError ){
985	<
986	<	sprintf( painCave.errMsg,
987	<	"There was an error setting the simulation information in fortran.\n" );
988	<	painCave.isFatal = 1;
989	<	painCave.severity = OPENMD_ERROR;
990	<	simError();
991	<	}
992	<
993	<
994	<	sprintf( checkPointMsg,
995	<	"succesfully sent the simulation information to fortran.\n");
996	<
997	<	errorCheckPoint();
998	<
999	<	// Setup number of neighbors in neighbor list if present
1000	<	if (simParams_->haveNeighborListNeighbors()) {
1001	<	int nlistNeighbors = simParams_->getNeighborListNeighbors();
1002	<	setNeighbors(&nlistNeighbors);
1003	<	}
1004	<
1005	<
1006	<	}
1007	<
1008	<
1009	<	void SimInfo::setupFortranParallel() {
1010	<	#ifdef IS_MPI
1011	<	//SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
1012	<	vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
1013	<	vector<int> localToGlobalCutoffGroupIndex;
1014	<	SimInfo::MoleculeIterator mi;
1015	<	Molecule::AtomIterator ai;
1016	<	Molecule::CutoffGroupIterator ci;
1017	<	Molecule* mol;
1018	<	Atom* atom;
1019	<	CutoffGroup* cg;
1020	<	mpiSimData parallelData;
1021	<	int isError;
1022	<
1023	<	for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
1024	<
1025	<	//local index(index in DataStorge) of atom is important
1026	<	for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
1027	<	localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
1028	<	}
1029	<
1030	<	//local index of cutoff group is trivial, it only depends on the order of travesing
1031	<	for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
1032	<	localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
1033	<	}
1034	<
1035	<	}
1036	<
1037	<	//fill up mpiSimData struct
1038	<	parallelData.nMolGlobal = getNGlobalMolecules();
1039	<	parallelData.nMolLocal = getNMolecules();
1040	<	parallelData.nAtomsGlobal = getNGlobalAtoms();
1041	<	parallelData.nAtomsLocal = getNAtoms();
1042	<	parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
1043	<	parallelData.nGroupsLocal = getNCutoffGroups();
1044	<	parallelData.myNode = worldRank;
1045	<	MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
1046	<
1047	<	//pass mpiSimData struct and index arrays to fortran
1048	<	setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
1049	<	&localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal),
1050	<	&localToGlobalCutoffGroupIndex[0], &isError);
1051	<
1052	<	if (isError) {
1053	<	sprintf(painCave.errMsg,
1054	<	"mpiRefresh errror: fortran didn't like something we gave it.\n");
1055	<	painCave.isFatal = 1;
1056	<	simError();
1057	<	}
1058	<
1059	<	sprintf(checkPointMsg, " mpiRefresh successful.\n");
1060	<	errorCheckPoint();
1061	<
1062	<	#endif
930	>	topologyDone_ = true;
931		}
932
1065	–
1066	–	void SimInfo::setupSwitchingFunction() {
1067	–
1068	–	}
1069	–
1070	–	void SimInfo::setupAccumulateBoxDipole() {
1071	–
1072	–	// we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
1073	–	if ( simParams_->haveAccumulateBoxDipole() )
1074	–	if ( simParams_->getAccumulateBoxDipole() ) {
1075	–	calcBoxDipole_ = true;
1076	–	}
1077	–
1078	–	}
1079	–
933		void SimInfo::addProperty(GenericData* genData) {
934		properties_.addProperty(genData);
935		}
#	Line 1111 \| Line 964 \| namespace OpenMD {
964		Molecule* mol;
965		RigidBody* rb;
966		Atom* atom;
967	+	CutoffGroup* cg;
968		SimInfo::MoleculeIterator mi;
969		Molecule::RigidBodyIterator rbIter;
970	<	Molecule::AtomIterator atomIter;;
970	>	Molecule::AtomIterator atomIter;
971	>	Molecule::CutoffGroupIterator cgIter;
972
973		for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
974
#	Line 1124 \| Line 979 \| namespace OpenMD {
979		for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
980		rb->setSnapshotManager(sman_);
981		}
982	+
983	+	for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) {
984	+	cg->setSnapshotManager(sman_);
985	+	}
986		}
987
988		}
#	Line 1352 \| Line 1211 \| namespace OpenMD {
1211
1212		det = intTensor.determinant();
1213		sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1214	<	volume = 4.0/3.0NumericConstant::PIpow(sysconstants,3.0/2.0)*sqrt(det);
1214	>	volume = 4.0/3.0NumericConstant::PIpow(sysconstants,geomCnst)*sqrt(det);
1215		return;
1216		}
1217
#	Line 1368 \| Line 1227 \| namespace OpenMD {
1227
1228		detI = intTensor.determinant();
1229		sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1230	<	volume = 4.0/3.0NumericConstant::PIpow(sysconstants,3.0/2.0)*sqrt(detI);
1230	>	volume = 4.0/3.0NumericConstant::PIpow(sysconstants,geomCnst)*sqrt(detI);
1231		return;
1232		}
1233		/*

Diff Legend

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

Comparing branches/development/src/brains/SimInfo.cpp (file contents): Revision 1532 by gezelter, Wed Dec 29 19:59:21 2010 UTC vs. Revision 1725 by gezelter, Sat May 26 18:13:43 2012 UTC

Diff Legend

Comparing branches/development/src/brains/SimInfo.cpp (file contents):
Revision 1532 by gezelter, Wed Dec 29 19:59:21 2010 UTC vs.
Revision 1725 by gezelter, Sat May 26 18:13:43 2012 UTC