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

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 749 by tim, Wed Nov 16 23:10:02 2005 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"
62	–	#include "UseTheForce/notifyCutoffs_interface.h"
64		#include "UseTheForce/DarkSide/switcheroo_interface.h"
65		#include "utils/MemoryUtils.hpp"
66		#include "utils/simError.h"
67		#include "selection/SelectionManager.hpp"
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 81 \| Line 85 \| namespace oopse {
85		return result;
86		}
87
88	<	SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
89	<	ForceField* ff, Globals* simParams) :
90	<	stamps_(stamps), forceField_(ff), simParams_(simParams),
87	<	ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
88	>	SimInfo::SimInfo(ForceField* ff, Globals* simParams) :
89	>	forceField_(ff), simParams_(simParams),
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	<
95	<	std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
98	>
99		MoleculeStamp* molStamp;
100		int nMolWithSameStamp;
101		int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
#	Line 100 \| Line 103 \| namespace oopse {
103		CutoffGroupStamp* cgStamp;
104		RigidBodyStamp* rbStamp;
105		int nRigidAtoms = 0;
106	<
107	<	for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
108	<	molStamp = i->first;
109	<	nMolWithSameStamp = i->second;
106	>
107	>	std::vector<Component*> components = simParams->getComponents();
108	>
109	>	for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
110	>	molStamp = (*i)->getMoleculeStamp();
111	>	nMolWithSameStamp = (*i)->getNMol();
112
113		addMoleculeStamp(molStamp, nMolWithSameStamp);
114
115		//calculate atoms in molecules
116		nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;
117
113	–
118		//calculate atoms in cutoff groups
119		int nAtomsInGroups = 0;
120		int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
121
122		for (int j=0; j < nCutoffGroupsInStamp; j++) {
123	<	cgStamp = molStamp->getCutoffGroup(j);
123	>	cgStamp = molStamp->getCutoffGroupStamp(j);
124		nAtomsInGroups += cgStamp->getNMembers();
125		}
126
#	Line 129 \| Line 133 \| namespace oopse {
133		int nRigidBodiesInStamp = molStamp->getNRigidBodies();
134
135		for (int j=0; j < nRigidBodiesInStamp; j++) {
136	<	rbStamp = molStamp->getRigidBody(j);
136	>	rbStamp = molStamp->getRigidBodyStamp(j);
137		nAtomsInRigidBodies += rbStamp->getNMembers();
138		}
139
#	Line 154 \| Line 158 \| namespace oopse {
158		+ nGlobalRigidBodies_;
159
160		nGlobalMols_ = molStampIds_.size();
157	–
158	–	#ifdef IS_MPI
161		molToProcMap_.resize(nGlobalMols_);
160	–	#endif
161	–
162		}
163
164		SimInfo::~SimInfo() {
#	Line 168 \| Line 168 \| namespace oopse {
168		}
169		molecules_.clear();
170
171	–	delete stamps_;
171		delete sman_;
172		delete simParams_;
173		delete forceField_;
#	Line 197 \| 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 222 \| 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 275 \| Line 276 \| namespace oopse {
276		}
277		}
278
279	<	}//end for (integrableObject)
280	<	}// end for (mol)
279	>	}
280	>	}
281
282		// n_constraints is local, so subtract them on each processor
283		ndf_local -= nConstraints_;
#	Line 293 \| Line 294 \| namespace oopse {
294
295		}
296
297	+	int SimInfo::getFdf() {
298	+	#ifdef IS_MPI
299	+	MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
300	+	#else
301	+	fdf_ = fdf_local;
302	+	#endif
303	+	return fdf_;
304	+	}
305	+
306		void SimInfo::calcNdfRaw() {
307		int ndfRaw_local;
308
#	Line 348 \| 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 413 \| Line 425 \| namespace oopse {
425		b = torsion->getAtomB()->getGlobalIndex();
426		c = torsion->getAtomC()->getGlobalIndex();
427		d = torsion->getAtomD()->getGlobalIndex();
428	+	std::set<int> rigidSetA = getRigidSet(a, atomGroups);
429	+	std::set<int> rigidSetB = getRigidSet(b, atomGroups);
430	+	std::set<int> rigidSetC = getRigidSet(c, atomGroups);
431	+	std::set<int> rigidSetD = getRigidSet(d, atomGroups);
432	+
433	+	exclude_.addPairs(rigidSetA, rigidSetB);
434	+	exclude_.addPairs(rigidSetA, rigidSetC);
435	+	exclude_.addPairs(rigidSetA, rigidSetD);
436	+	exclude_.addPairs(rigidSetB, rigidSetC);
437	+	exclude_.addPairs(rigidSetB, rigidSetD);
438	+	exclude_.addPairs(rigidSetC, rigidSetD);
439	+
440	+	/*
441	+	exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
442	+	exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
443	+	exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
444	+	exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
445	+	exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
446	+	exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
447	+
448	+
449	+	exclude_.addPair(a, b);
450	+	exclude_.addPair(a, c);
451	+	exclude_.addPair(a, d);
452	+	exclude_.addPair(b, c);
453	+	exclude_.addPair(b, d);
454	+	exclude_.addPair(c, d);
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);
#	Line 460 \| 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 525 \| Line 575 \| namespace oopse {
575		b = torsion->getAtomB()->getGlobalIndex();
576		c = torsion->getAtomC()->getGlobalIndex();
577		d = torsion->getAtomD()->getGlobalIndex();
578	+
579	+	std::set<int> rigidSetA = getRigidSet(a, atomGroups);
580	+	std::set<int> rigidSetB = getRigidSet(b, atomGroups);
581	+	std::set<int> rigidSetC = getRigidSet(c, atomGroups);
582	+	std::set<int> rigidSetD = getRigidSet(d, atomGroups);
583	+
584	+	exclude_.removePairs(rigidSetA, rigidSetB);
585	+	exclude_.removePairs(rigidSetA, rigidSetC);
586	+	exclude_.removePairs(rigidSetA, rigidSetD);
587	+	exclude_.removePairs(rigidSetB, rigidSetC);
588	+	exclude_.removePairs(rigidSetB, rigidSetD);
589	+	exclude_.removePairs(rigidSetC, rigidSetD);
590	+
591	+	/*
592	+	exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
593	+	exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
594	+	exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
595	+	exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
596	+	exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
597	+	exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
598	+
599	+
600	+	exclude_.removePair(a, b);
601	+	exclude_.removePair(a, c);
602	+	exclude_.removePair(a, d);
603	+	exclude_.removePair(b, c);
604	+	exclude_.removePair(b, d);
605	+	exclude_.removePair(c, d);
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);
#	Line 594 \| 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 603 \| Line 692 \| namespace oopse {
692		painCave.isFatal = 1;
693		simError();
694		}
606	–
607	–
608	–	setupCutoff();
695
696		calcNdf();
697		calcNdfRaw();
#	Line 655 \| 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 743 \| 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 762 \| 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	<
768	<	if (simParams_->haveDielectric()) {
769	<	fInfo_.dielect = simParams_->getDielectric();
770	<	} else {
771	<	sprintf(painCave.errMsg,
772	<	"SimSetup Error: No Dielectric constant was set.\n"
773	<	"\tYou are trying to use Reaction Field without"
774	<	"\tsetting a dielectric constant!\n");
775	<	painCave.isFatal = 1;
776	<	simError();
777	<	}
778	<	}
779	<
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 793 \| 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 843 \| 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
867	–	#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 916 \| Line 1017 \| namespace oopse {
1017		}
1018
1019		sprintf(checkPointMsg, " mpiRefresh successful.\n");
1020	<	MPIcheckPoint();
1020	>	errorCheckPoint();
1021
1022	<
1022	>	#endif
1023		}
1024
1025	<	#endif
1025	>	void SimInfo::setupCutoff() {
1026	>
1027	>	ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
1028
1029	<	double SimInfo::calcMaxCutoffRadius() {
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::set<AtomType*> atomTypes;
1037	<	std::set<AtomType*>::iterator i;
1038	<	std::vector<double> cutoffRadius;
1039	<
1040	<	//get the unique atom types
1041	<	atomTypes = getUniqueAtomTypes();
935	<
936	<	//query the max cutoff radius among these atom types
937	<	for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
938	<	cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
939	<	}
940	<
941	<	double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
942	<	#ifdef IS_MPI
943	<	//pick the max cutoff radius among the processors
944	<	#endif
945	<
946	<	return maxCutoffRadius;
947	<	}
948	<
949	<	void SimInfo::getCutoff(double& rcut, double& rsw) {
950	<
951	<	if (fInfo_.SIM_uses_Charges \| fInfo_.SIM_uses_Dipoles \| fInfo_.SIM_uses_RF) {
952	<
953	<	if (!simParams_->haveCutoffRadius()){
954	<	sprintf(painCave.errMsg,
955	<	"SimCreator Warning: No value was set for the cutoffRadius.\n"
956	<	"\tOOPSE will use a default value of 15.0 angstroms"
957	<	"\tfor the cutoffRadius.\n");
958	<	painCave.isFatal = 0;
959	<	simError();
960	<	rcut = 15.0;
961	<	} else{
962	<	rcut = simParams_->getCutoffRadius();
963	<	}
964	<
965	<	if (!simParams_->haveSwitchingRadius()){
966	<	sprintf(painCave.errMsg,
967	<	"SimCreator Warning: No value was set for switchingRadius.\n"
968	<	"\tOOPSE will use a default value of\n"
969	<	"\t0.85 * cutoffRadius for the switchingRadius\n");
970	<	painCave.isFatal = 0;
971	<	simError();
972	<	rsw = 0.85 * rcut;
973	<	} else{
974	<	rsw = simParams_->getSwitchingRadius();
975	<	}
976	<
977	<	} else {
978	<	// if charge, dipole or reaction field is not used and the cutofff radius is not specified in
979	<	//meta-data file, the maximum cutoff radius calculated from forcefiled will be used
980	<
981	<	if (simParams_->haveCutoffRadius()) {
982	<	rcut = simParams_->getCutoffRadius();
983	<	} else {
984	<	//set cutoff radius to the maximum cutoff radius based on atom types in the whole system
985	<	rcut = calcMaxCutoffRadius();
986	<	}
987	<
988	<	if (simParams_->haveSwitchingRadius()) {
989	<	rsw = simParams_->getSwitchingRadius();
990	<	} else {
991	<	rsw = rcut;
992	<	}
993	<
994	<	}
995	<	}
1036	>	std::string myPolicy;
1037	>	if (forceFieldOptions_.haveCutoffPolicy()){
1038	>	myPolicy = forceFieldOptions_.getCutoffPolicy();
1039	>	}else if (simParams_->haveCutoffPolicy()) {
1040	>	myPolicy = simParams_->getCutoffPolicy();
1041	>	}
1042
1043	<	void SimInfo::setupCutoff() {
998	<	getCutoff(rcut_, rsw_);
999	<	double rnblist = rcut_ + 1; // skin of neighbor list
1000	<
1001	<	//Pass these cutoff radius etc. to fortran. This function should be called once and only once
1002	<
1003	<	int cp = TRADITIONAL_CUTOFF_POLICY;
1004	<	if (simParams_->haveCutoffPolicy()) {
1005	<	std::string myPolicy = simParams_->getCutoffPolicy();
1043	>	if (!myPolicy.empty()){
1044		toUpper(myPolicy);
1045		if (myPolicy == "MIX") {
1046		cp = MIX_CUTOFF_POLICY;
#	Line 1021 \| Line 1059 \| namespace oopse {
1059		}
1060		}
1061		}
1062	<	}
1062	>	}
1063	>	notifyFortranCutoffPolicy(&cp);
1064
1065	<
1065	>	// Check the Skin Thickness for neighborlists
1066	>	RealType skin;
1067		if (simParams_->haveSkinThickness()) {
1068	<	double skinThickness = simParams_->getSkinThickness();
1069	<	}
1068	>	skin = simParams_->getSkinThickness();
1069	>	notifyFortranSkinThickness(&skin);
1070	>	}
1071	>
1072	>	// Check if the cutoff was set explicitly:
1073	>	if (simParams_->haveCutoffRadius()) {
1074	>	rcut_ = simParams_->getCutoffRadius();
1075	>	if (simParams_->haveSwitchingRadius()) {
1076	>	rsw_ = simParams_->getSwitchingRadius();
1077	>	} else {
1078	>	if (fInfo_.SIM_uses_Charges \|
1079	>	fInfo_.SIM_uses_Dipoles \|
1080	>	fInfo_.SIM_uses_RF) {
1081	>
1082	>	rsw_ = 0.85 * rcut_;
1083	>	sprintf(painCave.errMsg,
1084	>	"SimCreator Warning: No value was set for the switchingRadius.\n"
1085	>	"\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n"
1086	>	"\tswitchingRadius = %f. for this simulation\n", rsw_);
1087	>	painCave.isFatal = 0;
1088	>	simError();
1089	>	} else {
1090	>	rsw_ = rcut_;
1091	>	sprintf(painCave.errMsg,
1092	>	"SimCreator Warning: No value was set for the switchingRadius.\n"
1093	>	"\tOOPSE will use the same value as the cutoffRadius.\n"
1094	>	"\tswitchingRadius = %f. for this simulation\n", rsw_);
1095	>	painCave.isFatal = 0;
1096	>	simError();
1097	>	}
1098	>	}
1099
1100	<	notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
1101	<	// also send cutoff notification to electrostatics
1102	<	setElectrostaticCutoffRadius(&rcut_, &rsw_);
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	>
1112	>	} else {
1113	>
1114	>	// For electrostatic atoms, we'll assume a large safe value:
1115	>	if (fInfo_.SIM_uses_Charges \| fInfo_.SIM_uses_Dipoles \| fInfo_.SIM_uses_RF) {
1116	>	sprintf(painCave.errMsg,
1117	>	"SimCreator Warning: No value was set for the cutoffRadius.\n"
1118	>	"\tOOPSE will use a default value of 15.0 angstroms"
1119	>	"\tfor the cutoffRadius.\n");
1120	>	painCave.isFatal = 0;
1121	>	simError();
1122	>	rcut_ = 15.0;
1123	>
1124	>	if (simParams_->haveElectrostaticSummationMethod()) {
1125	>	std::string myMethod = simParams_->getElectrostaticSummationMethod();
1126	>	toUpper(myMethod);
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"
1139	>	"\teven though the electrostaticSummationMethod was\n"
1140	>	"\tset to %s\n", myMethod.c_str());
1141	>	painCave.isFatal = 1;
1142	>	simError();
1143	>	}
1144	>	}
1145	>	}
1146	>
1147	>	if (simParams_->haveSwitchingRadius()){
1148	>	rsw_ = simParams_->getSwitchingRadius();
1149	>	} else {
1150	>	sprintf(painCave.errMsg,
1151	>	"SimCreator Warning: No value was set for switchingRadius.\n"
1152	>	"\tOOPSE will use a default value of\n"
1153	>	"\t0.85 * cutoffRadius for the switchingRadius\n");
1154	>	painCave.isFatal = 0;
1155	>	simError();
1156	>	rsw_ = 0.85 * rcut_;
1157	>	}
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.
1164	>
1165	>	notifyFortranYouAreOnYourOwn();
1166	>
1167	>	}
1168	>	}
1169		}
1170
1171		void SimInfo::setupElectrostaticSummationMethod( int isError ) {
#	Line 1038 \| 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;
1045	–	alphaVal = simParams_->getDampingAlpha();
1046	–	dielectric = simParams_->getDielectric();
1180
1181		if (simParams_->haveElectrostaticSummationMethod()) {
1182		std::string myMethod = simParams_->getElectrostaticSummationMethod();
#	Line 1060 \| 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,
1210	<	"SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() );
1210	>	"SimInfo error: Unknown electrostaticSummationMethod.\n"
1211	>	"\t(Input file specified %s .)\n"
1212	>	"\telectrostaticSummationMethod must be one of: \"none\",\n"
1213	>	"\t\"shifted_potential\", \"shifted_force\", or \n"
1214	>	"\t\"reaction_field\".\n", myMethod.c_str() );
1215		painCave.isFatal = 1;
1216		simError();
1217		}
#	Line 1084 \| 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. A default value of %f (1/ang) will be used.", alphaVal);
1241	>	"SimInfo warning: dampingAlpha was not specified in the input file.\n"
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,
1252	<	"SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() );
1252	>	"SimInfo error: Unknown electrostaticScreeningMethod.\n"
1253	>	"\t(Input file specified %s .)\n"
1254	>	"\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
1255	>	"or \"damped\".\n", myScreen.c_str() );
1256		painCave.isFatal = 1;
1257		simError();
1258		}
#	Line 1102 \| Line 1261 \| namespace oopse {
1261
1262		// let's pass some summation method variables to fortran
1263		setElectrostaticSummationMethod( &esm );
1264	+	setFortranElectrostaticMethod( &esm );
1265		setScreeningMethod( &sm );
1266		setDampingAlpha( &alphaVal );
1267		setReactionFieldDielectric( &dielectric );
1268	<	initFortranFF( &esm, &errorOut );
1268	>	initFortranFF( &errorOut );
1269		}
1270
1271		void SimInfo::setupSwitchingFunction() {
#	Line 1134 \| 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 1190 \| 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 1216 \| 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 1252 \| 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 1287 \| 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 1304 \| 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 1342 \| 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 1364 \| 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 1377 \| 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 749 by tim, Wed Nov 16 23:10:02 2005 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 749 by tim, Wed Nov 16 23:10:02 2005 UTC vs.
Revision 1277 by gezelter, Mon Jul 14 12:35:58 2008 UTC