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

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 608 by chrisfen, Fri Sep 16 21:07:45 2005 UTC vs.
Revision 809 by gezelter, Mon Dec 12 19:32:50 2005 UTC

#	Line 48 \| Line 48
48
49		#include <algorithm>
50		#include <set>
51	+	#include <map>
52
53		#include "brains/SimInfo.hpp"
54		#include "math/Vector3.hpp"
55		#include "primitives/Molecule.hpp"
56		#include "UseTheForce/fCutoffPolicy.h"
57		#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
58	+	#include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
59	+	#include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
60		#include "UseTheForce/doForces_interface.h"
61	<	#include "UseTheForce/notifyCutoffs_interface.h"
61	>	#include "UseTheForce/DarkSide/electrostatic_interface.h"
62	>	#include "UseTheForce/DarkSide/switcheroo_interface.h"
63		#include "utils/MemoryUtils.hpp"
64		#include "utils/simError.h"
65		#include "selection/SelectionManager.hpp"
#	Line 66 \| Line 70 \| namespace oopse {
70		#endif
71
72		namespace oopse {
73	+	std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) {
74	+	std::map<int, std::set<int> >::iterator i = container.find(index);
75	+	std::set<int> result;
76	+	if (i != container.end()) {
77	+	result = i->second;
78	+	}
79
80	<	SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
81	<	ForceField* ff, Globals* simParams) :
82	<	stamps_(stamps), forceField_(ff), simParams_(simParams),
80	>	return result;
81	>	}
82	>
83	>	SimInfo::SimInfo(ForceField* ff, Globals* simParams) :
84	>	forceField_(ff), simParams_(simParams),
85		ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
86		nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
87		nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
#	Line 77 \| Line 89 \| namespace oopse {
89		nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0),
90		sman_(NULL), fortranInitialized_(false) {
91
80	–
81	–	std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
92		MoleculeStamp* molStamp;
93		int nMolWithSameStamp;
94		int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
95	<	int nGroups = 0; //total cutoff groups defined in meta-data file
95	>	int nGroups = 0; //total cutoff groups defined in meta-data file
96		CutoffGroupStamp* cgStamp;
97		RigidBodyStamp* rbStamp;
98		int nRigidAtoms = 0;
99	<
100	<	for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
101	<	molStamp = i->first;
102	<	nMolWithSameStamp = i->second;
99	>	std::vector<Component*> components = simParams->getComponents();
100	>
101	>	for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
102	>	molStamp = (*i)->getMoleculeStamp();
103	>	nMolWithSameStamp = (*i)->getNMol();
104
105		addMoleculeStamp(molStamp, nMolWithSameStamp);
106
107		//calculate atoms in molecules
108		nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;
109
99	–
110		//calculate atoms in cutoff groups
111		int nAtomsInGroups = 0;
112		int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
113
114		for (int j=0; j < nCutoffGroupsInStamp; j++) {
115	<	cgStamp = molStamp->getCutoffGroup(j);
115	>	cgStamp = molStamp->getCutoffGroupStamp(j);
116		nAtomsInGroups += cgStamp->getNMembers();
117		}
118
119		nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
120	+
121		nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;
122
123		//calculate atoms in rigid bodies
#	Line 114 \| Line 125 \| namespace oopse {
125		int nRigidBodiesInStamp = molStamp->getNRigidBodies();
126
127		for (int j=0; j < nRigidBodiesInStamp; j++) {
128	<	rbStamp = molStamp->getRigidBody(j);
128	>	rbStamp = molStamp->getRigidBodyStamp(j);
129		nAtomsInRigidBodies += rbStamp->getNMembers();
130		}
131
#	Line 123 \| Line 134 \| namespace oopse {
134
135		}
136
137	<	//every free atom (atom does not belong to cutoff groups) is a cutoff group
138	<	//therefore the total number of cutoff groups in the system is equal to
139	<	//the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
140	<	//file plus the number of cutoff groups defined in meta-data file
137	>	//every free atom (atom does not belong to cutoff groups) is a cutoff
138	>	//group therefore the total number of cutoff groups in the system is
139	>	//equal to the total number of atoms minus number of atoms belong to
140	>	//cutoff group defined in meta-data file plus the number of cutoff
141	>	//groups defined in meta-data file
142		nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
143
144	<	//every free atom (atom does not belong to rigid bodies) is an integrable object
145	<	//therefore the total number of integrable objects in the system is equal to
146	<	//the total number of atoms minus number of atoms belong to rigid body defined in meta-data
147	<	//file plus the number of rigid bodies defined in meta-data file
148	<	nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
149	<
144	>	//every free atom (atom does not belong to rigid bodies) is an
145	>	//integrable object therefore the total number of integrable objects
146	>	//in the system is equal to the total number of atoms minus number of
147	>	//atoms belong to rigid body defined in meta-data file plus the number
148	>	//of rigid bodies defined in meta-data file
149	>	nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
150	>	+ nGlobalRigidBodies_;
151	>
152		nGlobalMols_ = molStampIds_.size();
153
154		#ifdef IS_MPI
#	Line 150 \| Line 164 \| namespace oopse {
164		}
165		molecules_.clear();
166
153	–	delete stamps_;
167		delete sman_;
168		delete simParams_;
169		delete forceField_;
#	Line 257 \| Line 270 \| namespace oopse {
270		}
271		}
272
273	<	}//end for (integrableObject)
274	<	}// end for (mol)
273	>	}
274	>	}
275
276		// n_constraints is local, so subtract them on each processor
277		ndf_local -= nConstraints_;
#	Line 337 \| Line 350 \| namespace oopse {
350		int b;
351		int c;
352		int d;
353	+
354	+	std::map<int, std::set<int> > atomGroups;
355	+
356	+	Molecule::RigidBodyIterator rbIter;
357	+	RigidBody* rb;
358	+	Molecule::IntegrableObjectIterator ii;
359	+	StuntDouble* integrableObject;
360
361	+	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
362	+	integrableObject = mol->nextIntegrableObject(ii)) {
363	+
364	+	if (integrableObject->isRigidBody()) {
365	+	rb = static_cast<RigidBody*>(integrableObject);
366	+	std::vector<Atom*> atoms = rb->getAtoms();
367	+	std::set<int> rigidAtoms;
368	+	for (int i = 0; i < atoms.size(); ++i) {
369	+	rigidAtoms.insert(atoms[i]->getGlobalIndex());
370	+	}
371	+	for (int i = 0; i < atoms.size(); ++i) {
372	+	atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
373	+	}
374	+	} else {
375	+	std::set<int> oneAtomSet;
376	+	oneAtomSet.insert(integrableObject->getGlobalIndex());
377	+	atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));
378	+	}
379	+	}
380	+
381	+
382	+
383		for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
384		a = bond->getAtomA()->getGlobalIndex();
385		b = bond->getAtomB()->getGlobalIndex();
#	Line 348 \| Line 390 \| namespace oopse {
390		a = bend->getAtomA()->getGlobalIndex();
391		b = bend->getAtomB()->getGlobalIndex();
392		c = bend->getAtomC()->getGlobalIndex();
393	+	std::set<int> rigidSetA = getRigidSet(a, atomGroups);
394	+	std::set<int> rigidSetB = getRigidSet(b, atomGroups);
395	+	std::set<int> rigidSetC = getRigidSet(c, atomGroups);
396
397	<	exclude_.addPair(a, b);
398	<	exclude_.addPair(a, c);
399	<	exclude_.addPair(b, c);
397	>	exclude_.addPairs(rigidSetA, rigidSetB);
398	>	exclude_.addPairs(rigidSetA, rigidSetC);
399	>	exclude_.addPairs(rigidSetB, rigidSetC);
400	>
401	>	//exclude_.addPair(a, b);
402	>	//exclude_.addPair(a, c);
403	>	//exclude_.addPair(b, c);
404		}
405
406		for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
#	Line 359 \| Line 408 \| namespace oopse {
408		b = torsion->getAtomB()->getGlobalIndex();
409		c = torsion->getAtomC()->getGlobalIndex();
410		d = torsion->getAtomD()->getGlobalIndex();
411	+	std::set<int> rigidSetA = getRigidSet(a, atomGroups);
412	+	std::set<int> rigidSetB = getRigidSet(b, atomGroups);
413	+	std::set<int> rigidSetC = getRigidSet(c, atomGroups);
414	+	std::set<int> rigidSetD = getRigidSet(d, atomGroups);
415
416	+	exclude_.addPairs(rigidSetA, rigidSetB);
417	+	exclude_.addPairs(rigidSetA, rigidSetC);
418	+	exclude_.addPairs(rigidSetA, rigidSetD);
419	+	exclude_.addPairs(rigidSetB, rigidSetC);
420	+	exclude_.addPairs(rigidSetB, rigidSetD);
421	+	exclude_.addPairs(rigidSetC, rigidSetD);
422	+
423	+	/*
424	+	exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
425	+	exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
426	+	exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
427	+	exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
428	+	exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
429	+	exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
430	+
431	+
432		exclude_.addPair(a, b);
433		exclude_.addPair(a, c);
434		exclude_.addPair(a, d);
435		exclude_.addPair(b, c);
436		exclude_.addPair(b, d);
437		exclude_.addPair(c, d);
438	+	*/
439		}
440
371	–	Molecule::RigidBodyIterator rbIter;
372	–	RigidBody* rb;
441		for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
442		std::vector<Atom*> atoms = rb->getAtoms();
443		for (int i = 0; i < atoms.size() -1 ; ++i) {
#	Line 394 \| Line 462 \| namespace oopse {
462		int b;
463		int c;
464		int d;
465	+
466	+	std::map<int, std::set<int> > atomGroups;
467	+
468	+	Molecule::RigidBodyIterator rbIter;
469	+	RigidBody* rb;
470	+	Molecule::IntegrableObjectIterator ii;
471	+	StuntDouble* integrableObject;
472
473	+	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
474	+	integrableObject = mol->nextIntegrableObject(ii)) {
475	+
476	+	if (integrableObject->isRigidBody()) {
477	+	rb = static_cast<RigidBody*>(integrableObject);
478	+	std::vector<Atom*> atoms = rb->getAtoms();
479	+	std::set<int> rigidAtoms;
480	+	for (int i = 0; i < atoms.size(); ++i) {
481	+	rigidAtoms.insert(atoms[i]->getGlobalIndex());
482	+	}
483	+	for (int i = 0; i < atoms.size(); ++i) {
484	+	atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
485	+	}
486	+	} else {
487	+	std::set<int> oneAtomSet;
488	+	oneAtomSet.insert(integrableObject->getGlobalIndex());
489	+	atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));
490	+	}
491	+	}
492	+
493	+
494		for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
495		a = bond->getAtomA()->getGlobalIndex();
496		b = bond->getAtomB()->getGlobalIndex();
#	Line 406 \| Line 502 \| namespace oopse {
502		b = bend->getAtomB()->getGlobalIndex();
503		c = bend->getAtomC()->getGlobalIndex();
504
505	<	exclude_.removePair(a, b);
506	<	exclude_.removePair(a, c);
507	<	exclude_.removePair(b, c);
505	>	std::set<int> rigidSetA = getRigidSet(a, atomGroups);
506	>	std::set<int> rigidSetB = getRigidSet(b, atomGroups);
507	>	std::set<int> rigidSetC = getRigidSet(c, atomGroups);
508	>
509	>	exclude_.removePairs(rigidSetA, rigidSetB);
510	>	exclude_.removePairs(rigidSetA, rigidSetC);
511	>	exclude_.removePairs(rigidSetB, rigidSetC);
512	>
513	>	//exclude_.removePair(a, b);
514	>	//exclude_.removePair(a, c);
515	>	//exclude_.removePair(b, c);
516		}
517
518		for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
#	Line 417 \| Line 521 \| namespace oopse {
521		c = torsion->getAtomC()->getGlobalIndex();
522		d = torsion->getAtomD()->getGlobalIndex();
523
524	+	std::set<int> rigidSetA = getRigidSet(a, atomGroups);
525	+	std::set<int> rigidSetB = getRigidSet(b, atomGroups);
526	+	std::set<int> rigidSetC = getRigidSet(c, atomGroups);
527	+	std::set<int> rigidSetD = getRigidSet(d, atomGroups);
528	+
529	+	exclude_.removePairs(rigidSetA, rigidSetB);
530	+	exclude_.removePairs(rigidSetA, rigidSetC);
531	+	exclude_.removePairs(rigidSetA, rigidSetD);
532	+	exclude_.removePairs(rigidSetB, rigidSetC);
533	+	exclude_.removePairs(rigidSetB, rigidSetD);
534	+	exclude_.removePairs(rigidSetC, rigidSetD);
535	+
536	+	/*
537	+	exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
538	+	exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
539	+	exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
540	+	exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
541	+	exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
542	+	exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
543	+
544	+
545		exclude_.removePair(a, b);
546		exclude_.removePair(a, c);
547		exclude_.removePair(a, d);
548		exclude_.removePair(b, c);
549		exclude_.removePair(b, d);
550		exclude_.removePair(c, d);
551	+	*/
552		}
553
428	–	Molecule::RigidBodyIterator rbIter;
429	–	RigidBody* rb;
554		for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
555		std::vector<Atom*> atoms = rb->getAtoms();
556		for (int i = 0; i < atoms.size() -1 ; ++i) {
#	Line 466 \| Line 590 \| namespace oopse {
590		int isError = 0;
591
592		setupElectrostaticSummationMethod( isError );
593	+	setupSwitchingFunction();
594
595		if(isError){
596		sprintf( painCave.errMsg,
#	Line 510 \| Line 635 \| namespace oopse {
635		int useLennardJones = 0;
636		int useElectrostatic = 0;
637		int useEAM = 0;
638	+	int useSC = 0;
639		int useCharge = 0;
640		int useDirectional = 0;
641		int useDipole = 0;
#	Line 521 \| Line 647 \| namespace oopse {
647		int useDirectionalAtom = 0;
648		int useElectrostatics = 0;
649		//usePBC and useRF are from simParams
650	<	int usePBC = simParams_->getPBC();
650	>	int usePBC = simParams_->getUsePeriodicBoundaryConditions();
651	>	int useRF;
652	>	int useSF;
653	>	std::string myMethod;
654
655	+	// set the useRF logical
656	+	useRF = 0;
657	+	useSF = 0;
658	+
659	+
660	+	if (simParams_->haveElectrostaticSummationMethod()) {
661	+	std::string myMethod = simParams_->getElectrostaticSummationMethod();
662	+	toUpper(myMethod);
663	+	if (myMethod == "REACTION_FIELD") {
664	+	useRF=1;
665	+	} else {
666	+	if (myMethod == "SHIFTED_FORCE") {
667	+	useSF = 1;
668	+	}
669	+	}
670	+	}
671	+
672		//loop over all of the atom types
673		for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
674		useLennardJones \|= (*i)->isLennardJones();
675		useElectrostatic \|= (*i)->isElectrostatic();
676		useEAM \|= (*i)->isEAM();
677	+	useSC \|= (*i)->isSC();
678		useCharge \|= (*i)->isCharge();
679		useDirectional \|= (*i)->isDirectional();
680		useDipole \|= (*i)->isDipole();
#	Line 578 \| Line 725 \| namespace oopse {
725		temp = useEAM;
726		MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
727
728	+	temp = useSC;
729	+	MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
730	+
731		temp = useShape;
732		MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
733
734		temp = useFLARB;
735		MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
736
737	+	temp = useRF;
738	+	MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
739	+
740	+	temp = useSF;
741	+	MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
742	+
743		#endif
744
745		fInfo_.SIM_uses_PBC = usePBC;
#	Line 596 \| Line 752 \| namespace oopse {
752		fInfo_.SIM_uses_StickyPower = useStickyPower;
753		fInfo_.SIM_uses_GayBerne = useGayBerne;
754		fInfo_.SIM_uses_EAM = useEAM;
755	+	fInfo_.SIM_uses_SC = useSC;
756		fInfo_.SIM_uses_Shapes = useShape;
757		fInfo_.SIM_uses_FLARB = useFLARB;
758	+	fInfo_.SIM_uses_RF = useRF;
759	+	fInfo_.SIM_uses_SF = useSF;
760
761	<	if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
762	<
761	>	if( myMethod == "REACTION_FIELD") {
762	>
763		if (simParams_->haveDielectric()) {
764		fInfo_.dielect = simParams_->getDielectric();
765		} else {
#	Line 610 \| Line 769 \| namespace oopse {
769		"\tsetting a dielectric constant!\n");
770		painCave.isFatal = 1;
771		simError();
772	<	}
614	<
615	<	} else {
616	<	fInfo_.dielect = 0.0;
772	>	}
773		}
774
775		}
#	Line 649 \| Line 805 \| namespace oopse {
805
806		totalMass = cg->getMass();
807		for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
808	<	mfact.push_back(atom->getMass()/totalMass);
808	>	// Check for massless groups - set mfact to 1 if true
809	>	if (totalMass != 0)
810	>	mfact.push_back(atom->getMass()/totalMass);
811	>	else
812	>	mfact.push_back( 1.0 );
813		}
814
815		}
#	Line 758 \| Line 918 \| namespace oopse {
918
919		#endif
920
921	<	double SimInfo::calcMaxCutoffRadius() {
762	<
763	<
764	<	std::set<AtomType*> atomTypes;
765	<	std::set<AtomType*>::iterator i;
766	<	std::vector<double> cutoffRadius;
767	<
768	<	//get the unique atom types
769	<	atomTypes = getUniqueAtomTypes();
770	<
771	<	//query the max cutoff radius among these atom types
772	<	for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
773	<	cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
774	<	}
775	<
776	<	double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
777	<	#ifdef IS_MPI
778	<	//pick the max cutoff radius among the processors
779	<	#endif
780	<
781	<	return maxCutoffRadius;
782	<	}
783	<
784	<	void SimInfo::getCutoff(double& rcut, double& rsw) {
921	>	void SimInfo::setupCutoff() {
922
923	<	if (fInfo_.SIM_uses_Charges \| fInfo_.SIM_uses_Dipoles \| fInfo_.SIM_uses_RF) {
787	<
788	<	if (!simParams_->haveRcut()){
789	<	sprintf(painCave.errMsg,
790	<	"SimCreator Warning: No value was set for the cutoffRadius.\n"
791	<	"\tOOPSE will use a default value of 15.0 angstroms"
792	<	"\tfor the cutoffRadius.\n");
793	<	painCave.isFatal = 0;
794	<	simError();
795	<	rcut = 15.0;
796	<	} else{
797	<	rcut = simParams_->getRcut();
798	<	}
799	<
800	<	if (!simParams_->haveRsw()){
801	<	sprintf(painCave.errMsg,
802	<	"SimCreator Warning: No value was set for switchingRadius.\n"
803	<	"\tOOPSE will use a default value of\n"
804	<	"\t0.95 * cutoffRadius for the switchingRadius\n");
805	<	painCave.isFatal = 0;
806	<	simError();
807	<	rsw = 0.95 * rcut;
808	<	} else{
809	<	rsw = simParams_->getRsw();
810	<	}
811	<
812	<	} else {
813	<	// if charge, dipole or reaction field is not used and the cutofff radius is not specified in
814	<	//meta-data file, the maximum cutoff radius calculated from forcefiled will be used
815	<
816	<	if (simParams_->haveRcut()) {
817	<	rcut = simParams_->getRcut();
818	<	} else {
819	<	//set cutoff radius to the maximum cutoff radius based on atom types in the whole system
820	<	rcut = calcMaxCutoffRadius();
821	<	}
822	<
823	<	if (simParams_->haveRsw()) {
824	<	rsw = simParams_->getRsw();
825	<	} else {
826	<	rsw = rcut;
827	<	}
828	<
829	<	}
830	<	}
831	<
832	<	void SimInfo::setupCutoff() {
833	<	getCutoff(rcut_, rsw_);
834	<	double rnblist = rcut_ + 1; // skin of neighbor list
835	<
836	<	//Pass these cutoff radius etc. to fortran. This function should be called once and only once
837	<
923	>	// Check the cutoff policy
924		int cp = TRADITIONAL_CUTOFF_POLICY;
925		if (simParams_->haveCutoffPolicy()) {
926		std::string myPolicy = simParams_->getCutoffPolicy();
927	+	toUpper(myPolicy);
928		if (myPolicy == "MIX") {
929		cp = MIX_CUTOFF_POLICY;
930		} else {
#	Line 854 \| Line 941 \| namespace oopse {
941		simError();
942		}
943		}
944	+	}
945	+	}
946	+	notifyFortranCutoffPolicy(&cp);
947	+
948	+	// Check the Skin Thickness for neighborlists
949	+	double skin;
950	+	if (simParams_->haveSkinThickness()) {
951	+	skin = simParams_->getSkinThickness();
952	+	notifyFortranSkinThickness(&skin);
953	+	}
954	+
955	+	// Check if the cutoff was set explicitly:
956	+	if (simParams_->haveCutoffRadius()) {
957	+	rcut_ = simParams_->getCutoffRadius();
958	+	if (simParams_->haveSwitchingRadius()) {
959	+	rsw_ = simParams_->getSwitchingRadius();
960	+	} else {
961	+	rsw_ = rcut_;
962	+	}
963	+	notifyFortranCutoffs(&rcut_, &rsw_);
964	+
965	+	} else {
966	+
967	+	// For electrostatic atoms, we'll assume a large safe value:
968	+	if (fInfo_.SIM_uses_Charges \| fInfo_.SIM_uses_Dipoles \| fInfo_.SIM_uses_RF) {
969	+	sprintf(painCave.errMsg,
970	+	"SimCreator Warning: No value was set for the cutoffRadius.\n"
971	+	"\tOOPSE will use a default value of 15.0 angstroms"
972	+	"\tfor the cutoffRadius.\n");
973	+	painCave.isFatal = 0;
974	+	simError();
975	+	rcut_ = 15.0;
976	+
977	+	if (simParams_->haveElectrostaticSummationMethod()) {
978	+	std::string myMethod = simParams_->getElectrostaticSummationMethod();
979	+	toUpper(myMethod);
980	+	if (myMethod == "SHIFTED_POTENTIAL" \|\| myMethod == "SHIFTED_FORCE") {
981	+	if (simParams_->haveSwitchingRadius()){
982	+	sprintf(painCave.errMsg,
983	+	"SimInfo Warning: A value was set for the switchingRadius\n"
984	+	"\teven though the electrostaticSummationMethod was\n"
985	+	"\tset to %s\n", myMethod.c_str());
986	+	painCave.isFatal = 1;
987	+	simError();
988	+	}
989	+	}
990	+	}
991	+
992	+	if (simParams_->haveSwitchingRadius()){
993	+	rsw_ = simParams_->getSwitchingRadius();
994	+	} else {
995	+	sprintf(painCave.errMsg,
996	+	"SimCreator Warning: No value was set for switchingRadius.\n"
997	+	"\tOOPSE will use a default value of\n"
998	+	"\t0.85 * cutoffRadius for the switchingRadius\n");
999	+	painCave.isFatal = 0;
1000	+	simError();
1001	+	rsw_ = 0.85 * rcut_;
1002	+	}
1003	+	notifyFortranCutoffs(&rcut_, &rsw_);
1004	+	} else {
1005	+	// We didn't set rcut explicitly, and we don't have electrostatic atoms, so
1006	+	// We'll punt and let fortran figure out the cutoffs later.
1007	+
1008	+	notifyFortranYouAreOnYourOwn();
1009	+
1010		}
1011		}
859	–	notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
1012		}
1013
1014		void SimInfo::setupElectrostaticSummationMethod( int isError ) {
1015
1016		int errorOut;
1017		int esm = NONE;
1018	+	int sm = UNDAMPED;
1019		double alphaVal;
1020	+	double dielectric;
1021
1022		errorOut = isError;
1023	+	alphaVal = simParams_->getDampingAlpha();
1024	+	dielectric = simParams_->getDielectric();
1025
1026		if (simParams_->haveElectrostaticSummationMethod()) {
1027		std::string myMethod = simParams_->getElectrostaticSummationMethod();
1028	+	toUpper(myMethod);
1029		if (myMethod == "NONE") {
1030		esm = NONE;
1031		} else {
1032	<	if (myMethod == "UNDAMPED_WOLF") {
1033	<	esm = UNDAMPED_WOLF;
1032	>	if (myMethod == "SWITCHING_FUNCTION") {
1033	>	esm = SWITCHING_FUNCTION;
1034		} else {
1035	<	if (myMethod == "DAMPED_WOLF") {
1036	<	esm = DAMPED_WOLF;
1037	<	if (!simParams_->haveDampingAlpha()) {
1038	<	//throw error
1039	<	sprintf( painCave.errMsg,
883	<	"SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", simParams_->getDampingAlpha());
884	<	painCave.isFatal = 0;
885	<	simError();
886	<	}
887	<	alphaVal = simParams_->getDampingAlpha();
888	<	} else {
889	<	if (myMethod == "REACTION_FIELD") {
890	<	esm = REACTION_FIELD;
1035	>	if (myMethod == "SHIFTED_POTENTIAL") {
1036	>	esm = SHIFTED_POTENTIAL;
1037	>	} else {
1038	>	if (myMethod == "SHIFTED_FORCE") {
1039	>	esm = SHIFTED_FORCE;
1040		} else {
1041	<	// throw error
1042	<	sprintf( painCave.errMsg,
1043	<	"SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() );
1044	<	painCave.isFatal = 1;
1045	<	simError();
1046	<	}
1047	<	}
1041	>	if (myMethod == "REACTION_FIELD") {
1042	>	esm = REACTION_FIELD;
1043	>	} else {
1044	>	// throw error
1045	>	sprintf( painCave.errMsg,
1046	>	"SimInfo error: Unknown electrostaticSummationMethod.\n"
1047	>	"\t(Input file specified %s .)\n"
1048	>	"\telectrostaticSummationMethod must be one of: \"none\",\n"
1049	>	"\t\"shifted_potential\", \"shifted_force\", or \n"
1050	>	"\t\"reaction_field\".\n", myMethod.c_str() );
1051	>	painCave.isFatal = 1;
1052	>	simError();
1053	>	}
1054	>	}
1055	>	}
1056		}
1057		}
1058		}
1059	<	initFortranFF( &esm, &alphaVal, &errorOut );
1059	>
1060	>	if (simParams_->haveElectrostaticScreeningMethod()) {
1061	>	std::string myScreen = simParams_->getElectrostaticScreeningMethod();
1062	>	toUpper(myScreen);
1063	>	if (myScreen == "UNDAMPED") {
1064	>	sm = UNDAMPED;
1065	>	} else {
1066	>	if (myScreen == "DAMPED") {
1067	>	sm = DAMPED;
1068	>	if (!simParams_->haveDampingAlpha()) {
1069	>	//throw error
1070	>	sprintf( painCave.errMsg,
1071	>	"SimInfo warning: dampingAlpha was not specified in the input file.\n"
1072	>	"\tA default value of %f (1/ang) will be used.\n", alphaVal);
1073	>	painCave.isFatal = 0;
1074	>	simError();
1075	>	}
1076	>	} else {
1077	>	// throw error
1078	>	sprintf( painCave.errMsg,
1079	>	"SimInfo error: Unknown electrostaticScreeningMethod.\n"
1080	>	"\t(Input file specified %s .)\n"
1081	>	"\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
1082	>	"or \"damped\".\n", myScreen.c_str() );
1083	>	painCave.isFatal = 1;
1084	>	simError();
1085	>	}
1086	>	}
1087	>	}
1088	>
1089	>	// let's pass some summation method variables to fortran
1090	>	setElectrostaticSumMethod( &esm );
1091	>	setFortranElectrostaticMethod( &esm );
1092	>	setScreeningMethod( &sm );
1093	>	setDampingAlpha( &alphaVal );
1094	>	setReactionFieldDielectric( &dielectric );
1095	>	initFortranFF( &errorOut );
1096		}
1097
1098	+	void SimInfo::setupSwitchingFunction() {
1099	+	int ft = CUBIC;
1100	+
1101	+	if (simParams_->haveSwitchingFunctionType()) {
1102	+	std::string funcType = simParams_->getSwitchingFunctionType();
1103	+	toUpper(funcType);
1104	+	if (funcType == "CUBIC") {
1105	+	ft = CUBIC;
1106	+	} else {
1107	+	if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1108	+	ft = FIFTH_ORDER_POLY;
1109	+	} else {
1110	+	// throw error
1111	+	sprintf( painCave.errMsg,
1112	+	"SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1113	+	painCave.isFatal = 1;
1114	+	simError();
1115	+	}
1116	+	}
1117	+	}
1118	+
1119	+	// send switching function notification to switcheroo
1120	+	setFunctionType(&ft);
1121	+
1122	+	}
1123	+
1124		void SimInfo::addProperty(GenericData* genData) {
1125		properties_.addProperty(genData);
1126		}

Diff Legend

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

Comparing trunk/src/brains/SimInfo.cpp (file contents): Revision 608 by chrisfen, Fri Sep 16 21:07:45 2005 UTC vs. Revision 809 by gezelter, Mon Dec 12 19:32:50 2005 UTC

Diff Legend

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 608 by chrisfen, Fri Sep 16 21:07:45 2005 UTC vs.
Revision 809 by gezelter, Mon Dec 12 19:32:50 2005 UTC