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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 407 by mmeineke, Wed Mar 26 20:22:02 2003 UTC vs.
Revision 489 by mmeineke, Thu Apr 10 20:08:56 2003 UTC

#	Line 66 \| Line 66 \| void SimSetup::createSim( void ){
66
67		MakeStamps *the_stamps;
68		Globals* the_globals;
69	+	ExtendedSystem* the_extendedsystem;
70		int i, j;
71
72		// get the stamps and globals;
#	Line 80 \| Line 81 \| void SimSetup::createSim( void ){
81		// get the ones we know are there, yet still may need some work.
82		n_components = the_globals->getNComponents();
83		strcpy( force_field, the_globals->getForceField() );
84	+
85	+	// get the ensemble and set up an extended system if we need it:
86		strcpy( ensemble, the_globals->getEnsemble() );
87	+	if( !strcasecmp( ensemble, "NPT" ) ) {
88	+	the_extendedsystem = new ExtendedSystem( simnfo );
89	+	the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
90	+	if (the_globals->haveTargetPressure())
91	+	the_extendedsystem->setTargetPressure(the_globals->getTargetPressure());
92	+	else {
93	+	sprintf( painCave.errMsg,
94	+	"SimSetup error: If you use the constant pressure\n"
95	+	" ensemble, you must set targetPressure.\n"
96	+	" This was found in the BASS file.\n");
97	+	painCave.isFatal = 1;
98	+	simError();
99	+	}
100	+
101	+	if (the_globals->haveTauThermostat())
102	+	the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
103	+	else if (the_globals->haveQmass())
104	+	the_extendedsystem->setQmass(the_globals->getQmass());
105	+	else {
106	+	sprintf( painCave.errMsg,
107	+	"SimSetup error: If you use one of the constant temperature\n"
108	+	" ensembles, you must set either tauThermostat or qMass.\n"
109	+	" Neither of these was found in the BASS file.\n");
110	+	painCave.isFatal = 1;
111	+	simError();
112	+	}
113	+
114	+	if (the_globals->haveTauBarostat())
115	+	the_extendedsystem->setTauBarostat(the_globals->getTauBarostat());
116	+	else {
117	+	sprintf( painCave.errMsg,
118	+	"SimSetup error: If you use the constant pressure\n"
119	+	" ensemble, you must set tauBarostat.\n"
120	+	" This was found in the BASS file.\n");
121	+	painCave.isFatal = 1;
122	+	simError();
123	+	}
124	+
125	+	} else if ( !strcasecmp( ensemble, "NVT") ) {
126	+	the_extendedsystem = new ExtendedSystem( simnfo );
127	+	the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
128	+
129	+	if (the_globals->haveTauThermostat())
130	+	the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
131	+	else if (the_globals->haveQmass())
132	+	the_extendedsystem->setQmass(the_globals->getQmass());
133	+	else {
134	+	sprintf( painCave.errMsg,
135	+	"SimSetup error: If you use one of the constant temperature\n"
136	+	" ensembles, you must set either tauThermostat or qMass.\n"
137	+	" Neither of these was found in the BASS file.\n");
138	+	painCave.isFatal = 1;
139	+	simError();
140	+	}
141	+
142	+	} else if ( !strcasecmp( ensemble, "NVE") ) {
143	+	} else {
144	+	sprintf( painCave.errMsg,
145	+	"SimSetup Warning. Unrecognized Ensemble -> %s, "
146	+	"reverting to NVE for this simulation.\n",
147	+	ensemble );
148	+	painCave.isFatal = 0;
149	+	simError();
150	+	strcpy( ensemble, "NVE" );
151	+	}
152		strcpy( simnfo->ensemble, ensemble );
153
154		strcpy( simnfo->mixingRule, the_globals->getMixingRule() );
155		simnfo->usePBC = the_globals->getPBC();
156
157	<
158	<
159	<	if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
160	<	else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
161	<	else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
162	<	else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
157	>	int usesDipoles = 0;
158	>	if( !strcmp( force_field, "TraPPE_Ex" ) ){
159	>	the_ff = new TraPPE_ExFF();
160	>	usesDipoles = 1;
161	>	}
162	>	else if( !strcasecmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
163		else{
164		sprintf( painCave.errMsg,
165		"SimSetup Error. Unrecognized force field -> %s\n",
#	Line 242 \| Line 310 \| void SimSetup::createSim( void ){
310
311		globalIndex = mpiSim->divideLabor();
312
245	–
246	–
313		// set up the local variables
314
315		int localMol, allMol;
316		int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
317	+
318	+	int* mol2proc = mpiSim->getMolToProcMap();
319	+	int* molCompType = mpiSim->getMolComponentType();
320
321		allMol = 0;
322		localMol = 0;
#	Line 259 \| Line 328 \| void SimSetup::createSim( void ){
328
329		for( j=0; j<components_nmol[i]; j++ ){
330
331	<	if( mpiSim->getMyMolStart() <= allMol &&
263	<	allMol <= mpiSim->getMyMolEnd() ){
331	>	if( mol2proc[allMol] == worldRank ){
332
333		local_atoms += comp_stamps[i]->getNAtoms();
334		local_bonds += comp_stamps[i]->getNBonds();
#	Line 279 \| Line 347 \| void SimSetup::createSim( void ){
347		if( local_atoms != simnfo->n_atoms ){
348		sprintf( painCave.errMsg,
349		"SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
350	<	" localAtom (%d) are note equal.\n",
350	>	" localAtom (%d) are not equal.\n",
351		simnfo->n_atoms,
352		local_atoms );
353		painCave.isFatal = 1;
#	Line 304 \| Line 372 \| void SimSetup::createSim( void ){
372		Atom::createArrays(simnfo->n_atoms);
373		the_atoms = new Atom*[simnfo->n_atoms];
374		the_molecules = new Molecule[simnfo->n_mol];
375	+	int molIndex;
376
377	+	// initialize the molecule's stampID's
378
379	+	#ifdef IS_MPI
380	+
381	+
382	+	molIndex = 0;
383	+	for(i=0; i<mpiSim->getTotNmol(); i++){
384	+
385	+	if(mol2proc[i] == worldRank ){
386	+	the_molecules[molIndex].setStampID( molCompType[i] );
387	+	the_molecules[molIndex].setMyIndex( molIndex );
388	+	the_molecules[molIndex].setGlobalIndex( i );
389	+	molIndex++;
390	+	}
391	+	}
392	+
393	+	#else // is_mpi
394	+
395	+	molIndex = 0;
396	+	for(i=0; i<n_components; i++){
397	+	for(j=0; j<components_nmol[i]; j++ ){
398	+	the_molecules[molIndex].setStampID( i );
399	+	the_molecules[molIndex].setMyIndex( molIndex );
400	+	the_molecules[molIndex].setGlobalIndex( molIndex );
401	+	molIndex++;
402	+	}
403	+	}
404	+
405	+
406	+	#endif // is_mpi
407	+
408	+
409		if( simnfo->n_SRI ){
410	<	the_sris = new SRI*[simnfo->n_SRI];
411	<	the_excludes = new int[2 * simnfo->n_SRI];
410	>
411	>	Exclude::createArray(simnfo->n_SRI);
412	>	the_excludes = new Exclude*[simnfo->n_SRI];
413	>	for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
414		simnfo->globalExcludes = new int;
415	<	simnfo->n_exclude = tot_SRI;
415	>	simnfo->n_exclude = simnfo->n_SRI;
416		}
417		else{
418
419	<	the_excludes = new int[2];
420	<	the_excludes[0] = 0;
421	<	the_excludes[1] = 0;
419	>	Exclude::createArray( 1 );
420	>	the_excludes = new Exclude*;
421	>	the_excludes[0] = new Exclude(0);
422	>	the_excludes[0]->setPair( 0,0 );
423		simnfo->globalExcludes = new int;
424		simnfo->globalExcludes[0] = 0;
425	<
323	<	simnfo->n_exclude = 1;
425	>	simnfo->n_exclude = 0;
426		}
427
428		// set the arrays into the SimInfo object
429
430		simnfo->atoms = the_atoms;
431	<	simnfo->sr_interactions = the_sris;
431	>	simnfo->molecules = the_molecules;
432		simnfo->nGlobalExcludes = 0;
433		simnfo->excludes = the_excludes;
434
#	Line 383 \| Line 485 \| void SimSetup::createSim( void ){
485
486		the_ff->setSimInfo( simnfo );
487
488	<	makeAtoms();
488	>	makeMolecules();
489		simnfo->identArray = new int[simnfo->n_atoms];
490	+	simnfo->molMembershipArray = new int[simnfo->n_atoms];
491		for(i=0; i<simnfo->n_atoms; i++){
492		simnfo->identArray[i] = the_atoms[i]->getIdent();
493		}
494	<
495	<	if( tot_bonds ){
393	<	makeBonds();
494	>	for(i=0; i< simnfo->n_mol; i++) {
495	>	the_molecules[i].atomicRollCall(simnfo->molMembershipArray);
496		}
497	<
396	<	if( tot_bends ){
397	<	makeBends();
398	<	}
399	<
400	<	if( tot_torsions ){
401	<	makeTorsions();
402	<	}
403	<
404	<
497	>
498		if (the_globals->getUseRF() ) {
499		simnfo->useReactionField = 1;
500
#	Line 443 \| Line 536 \| void SimSetup::createSim( void ){
536		}
537		simnfo->dielectric = the_globals->getDielectric();
538		} else {
539	<	if (simnfo->n_dipoles) {
539	>	if (usesDipoles) {
540
541		if( !the_globals->haveECR() ){
542		sprintf( painCave.errMsg,
543	<	"SimSetup Warning: using default value of 1/2 the smallest"
543	>	"SimSetup Warning: using default value of 1/2 the smallest "
544		"box length for the electrostaticCutoffRadius.\n"
545		"I hope you have a very fast processor!\n");
546		painCave.isFatal = 0;
#	Line 463 \| Line 556 \| void SimSetup::createSim( void ){
556
557		if( !the_globals->haveEST() ){
558		sprintf( painCave.errMsg,
559	<	"SimSetup Warning: using default value of 5% of the"
559	>	"SimSetup Warning: using default value of 5%% of the "
560		"electrostaticCutoffRadius for the "
561		"electrostaticSkinThickness\n"
562		);
#	Line 638 \| Line 731 \| void SimSetup::createSim( void ){
731
732		// new AllLong( simnfo );
733
641	–	if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
642	–	if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff );
643	–	if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff );
644	–	if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
734
735	<
735	>	if( !strcmp( force_field, "TraPPE_Ex" ) ){
736	>	new Symplectic(simnfo, the_ff, the_extendedsystem);
737	>	}
738	>	else if( !strcmp( force_field, "LJ" ) ){
739	>	new Verlet( *simnfo, the_ff, the_extendedsystem );
740	>	}
741	>	else {
742	>	std::cerr << "I'm a bug.\n";
743	>	fprintf( stderr, "Ima bug. stderr %s\n", force_field);
744	>	}
745	>	#ifdef IS_MPI
746	>	mpiSim->mpiRefresh();
747	>	#endif
748
749		// initialize the Fortran
750	<
750	>
751	>
752		simnfo->refreshSim();
753
754		if( !strcmp( simnfo->mixingRule, "standard") ){
#	Line 674 \| Line 776 \| void SimSetup::makeMolecules( void ){
776
777		void SimSetup::makeMolecules( void ){
778
779	<	int i, j, exI, exJ, tempEx, stampID, atomOffset;
779	>	int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
780		molInit info;
781		DirectionalAtom* dAtom;
782	+	LinkedAssign* extras;
783	+	LinkedAssign* current_extra;
784		AtomStamp* currentAtom;
785		BondStamp* currentBond;
786		BendStamp* currentBend;
787		TorsionStamp* currentTorsion;
788	+
789	+	bond_pair* theBonds;
790	+	bend_set* theBends;
791	+	torsion_set* theTorsions;
792	+
793
794		//init the forceField paramters
795
796		the_ff->readParams();
797
798
799	<	// init the molecules
799	>	// init the atoms
800
801	+	double ux, uy, uz, u, uSqr;
802	+
803		atomOffset = 0;
804	+	excludeOffset = 0;
805		for(i=0; i<simnfo->n_mol; i++){
806
807		stampID = the_molecules[i].getStampID();
#	Line 698 \| Line 810 \| void SimSetup::makeMolecules( void ){
810		info.nBonds = comp_stamps[stampID]->getNBonds();
811		info.nBends = comp_stamps[stampID]->getNBends();
812		info.nTorsions = comp_stamps[stampID]->getNTorsions();
813	<
813	>	info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
814	>
815		info.myAtoms = &the_atoms[atomOffset];
816	+	info.myExcludes = &the_excludes[excludeOffset];
817		info.myBonds = new Bond*[info.nBonds];
818		info.myBends = new Bend*[info.nBends];
819	<	info.myTorsions = new Torsions*[info.nTorsions];
819	>	info.myTorsions = new Torsion*[info.nTorsions];
820
821		theBonds = new bond_pair[info.nBonds];
822		theBends = new bend_set[info.nBends];
#	Line 712 \| Line 826 \| void SimSetup::makeMolecules( void ){
826
827		for(j=0; j<info.nAtoms; j++){
828
829	<	currentAtom = theComponents[stampID]->getAtom( j );
829	>	currentAtom = comp_stamps[stampID]->getAtom( j );
830		if( currentAtom->haveOrientation() ){
831
832		dAtom = new DirectionalAtom(j + atomOffset);
#	Line 747 \| Line 861 \| void SimSetup::makeMolecules( void ){
861		}
862
863		// make the bonds
864	<	for(j=0; j<nBonds; j++){
864	>	for(j=0; j<info.nBonds; j++){
865
866		currentBond = comp_stamps[stampID]->getBond( j );
867		theBonds[j].a = currentBond->getA() + atomOffset;
868		theBonds[j].b = currentBond->getB() + atomOffset;
869
870	<	exI = theBonds[i].a;
871	<	exJ = theBonds[i].b;
870	>	exI = theBonds[j].a;
871	>	exJ = theBonds[j].b;
872
873		// exclude_I must always be the smaller of the pair
874		if( exI > exJ ){
#	Line 763 \| Line 877 \| void SimSetup::makeMolecules( void ){
877		exJ = tempEx;
878		}
879		#ifdef IS_MPI
880	+	tempEx = exI;
881	+	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
882	+	tempEx = exJ;
883	+	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
884
885	<	the_excludes[index*2] =
768	<	the_atoms[exI]->getGlobalIndex() + 1;
769	<	the_excludes[index*2 + 1] =
770	<	the_atoms[exJ]->getGlobalIndex() + 1;
771	<
885	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
886		#else // isn't MPI
773	–
774	–	the_excludes[index*2] = exI + 1;
775	–	the_excludes[index*2 + 1] = exJ + 1;
776	–	// fortran index from 1 (hence the +1 in the indexing)
887
888	+	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
889		#endif //is_mpi
890	<
891	<	}
890	>	}
891	>	excludeOffset += info.nBonds;
892
893	<
894	<
895	<
896	<
897	<
898	<
899	<
789	<
790	<
791	<
792	<
793	<
794	<
795	<
796	<	void SimSetup::makeAtoms( void ){
797	<
798	<	int i, j, k, index;
799	<	double ux, uy, uz, uSqr, u;
800	<	AtomStamp* current_atom;
801	<
802	<	DirectionalAtom* dAtom;
803	<	int molIndex, molStart, molEnd, nMemb, lMolIndex;
804	<
805	<	lMolIndex = 0;
806	<	molIndex = 0;
807	<	index = 0;
808	<	for( i=0; i<n_components; i++ ){
809	<
810	<	for( j=0; j<components_nmol[i]; j++ ){
811	<
812	<	#ifdef IS_MPI
813	<	if( mpiSim->getMyMolStart() <= molIndex &&
814	<	molIndex <= mpiSim->getMyMolEnd() ){
815	<	#endif // is_mpi
816	<
817	<	molStart = index;
818	<	nMemb = comp_stamps[i]->getNAtoms();
819	<	for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
893	>	//make the bends
894	>	for(j=0; j<info.nBends; j++){
895	>
896	>	currentBend = comp_stamps[stampID]->getBend( j );
897	>	theBends[j].a = currentBend->getA() + atomOffset;
898	>	theBends[j].b = currentBend->getB() + atomOffset;
899	>	theBends[j].c = currentBend->getC() + atomOffset;
900
901	<	current_atom = comp_stamps[i]->getAtom( k );
822	<	if( current_atom->haveOrientation() ){
901	>	if( currentBend->haveExtras() ){
902
903	<	dAtom = new DirectionalAtom(index);
904	<	simnfo->n_oriented++;
826	<	the_atoms[index] = dAtom;
903	>	extras = currentBend->getExtras();
904	>	current_extra = extras;
905
906	<	ux = current_atom->getOrntX();
907	<	uy = current_atom->getOrntY();
908	<	uz = current_atom->getOrntZ();
909	<
910	<	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
911	<
912	<	u = sqrt( uSqr );
913	<	ux = ux / u;
914	<	uy = uy / u;
915	<	uz = uz / u;
916	<
917	<	dAtom->setSUx( ux );
918	<	dAtom->setSUy( uy );
919	<	dAtom->setSUz( uz );
920	<	}
921	<	else{
922	<	the_atoms[index] = new GeneralAtom(index);
923	<	}
924	<	the_atoms[index]->setType( current_atom->getType() );
925	<	the_atoms[index]->setIndex( index );
926	<
927	<	// increment the index and repeat;
928	<	index++;
929	<	}
930	<
931	<	molEnd = index -1;
932	<	the_molecules[lMolIndex].setNMembers( nMemb );
855	<	the_molecules[lMolIndex].setStartAtom( molStart );
856	<	the_molecules[lMolIndex].setEndAtom( molEnd );
857	<	the_molecules[lMolIndex].setStampID( i );
858	<	lMolIndex++;
859	<
860	<	#ifdef IS_MPI
861	<	}
862	<	#endif //is_mpi
863	<
864	<	molIndex++;
865	<	}
866	<	}
867	<
868	<	#ifdef IS_MPI
869	<	for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
870	<
871	<	delete[] globalIndex;
872	<
873	<	mpiSim->mpiRefresh();
874	<	#endif //IS_MPI
875	<
876	<	the_ff->initializeAtoms();
877	<	}
878	<
879	<	void SimSetup::makeBonds( void ){
880	<
881	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
882	<	bond_pair* the_bonds;
883	<	BondStamp* current_bond;
884	<
885	<	the_bonds = new bond_pair[tot_bonds];
886	<	index = 0;
887	<	offset = 0;
888	<	molIndex = 0;
889	<
890	<	for( i=0; i<n_components; i++ ){
891	<
892	<	for( j=0; j<components_nmol[i]; j++ ){
893	<
894	<	#ifdef IS_MPI
895	<	if( mpiSim->getMyMolStart() <= molIndex &&
896	<	molIndex <= mpiSim->getMyMolEnd() ){
897	<	#endif // is_mpi
898	<
899	<	for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
906	>	while( current_extra != NULL ){
907	>	if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
908	>
909	>	switch( current_extra->getType() ){
910	>
911	>	case 0:
912	>	theBends[j].ghost =
913	>	current_extra->getInt() + atomOffset;
914	>	theBends[j].isGhost = 1;
915	>	break;
916	>
917	>	case 1:
918	>	theBends[j].ghost =
919	>	(int)current_extra->getDouble() + atomOffset;
920	>	theBends[j].isGhost = 1;
921	>	break;
922	>
923	>	default:
924	>	sprintf( painCave.errMsg,
925	>	"SimSetup Error: ghostVectorSource was neither a "
926	>	"double nor an int.\n"
927	>	"-->Bend[%d] in %s\n",
928	>	j, comp_stamps[stampID]->getID() );
929	>	painCave.isFatal = 1;
930	>	simError();
931	>	}
932	>	}
933
934	<	current_bond = comp_stamps[i]->getBond( k );
935	<	the_bonds[index].a = current_bond->getA() + offset;
936	<	the_bonds[index].b = current_bond->getB() + offset;
937	<
938	<	exI = the_bonds[index].a;
939	<	exJ = the_bonds[index].b;
940	<
941	<	// exclude_I must always be the smaller of the pair
942	<	if( exI > exJ ){
910	<	tempEx = exI;
911	<	exI = exJ;
912	<	exJ = tempEx;
934	>	else{
935	>
936	>	sprintf( painCave.errMsg,
937	>	"SimSetup Error: unhandled bend assignment:\n"
938	>	" -->%s in Bend[%d] in %s\n",
939	>	current_extra->getlhs(),
940	>	j, comp_stamps[stampID]->getID() );
941	>	painCave.isFatal = 1;
942	>	simError();
943		}
914	–
944
945	<	#ifdef IS_MPI
917	<
918	<	the_excludes[index*2] =
919	<	the_atoms[exI]->getGlobalIndex() + 1;
920	<	the_excludes[index*2 + 1] =
921	<	the_atoms[exJ]->getGlobalIndex() + 1;
922	<
923	<	#else // isn't MPI
924	<
925	<	the_excludes[index*2] = exI + 1;
926	<	the_excludes[index*2 + 1] = exJ + 1;
927	<	// fortran index from 1 (hence the +1 in the indexing)
928	<	#endif //is_mpi
929	<
930	<	// increment the index and repeat;
931	<	index++;
945	>	current_extra = current_extra->getNext();
946		}
947	<	offset += comp_stamps[i]->getNAtoms();
947	>	}
948	>
949	>	if( !theBends[j].isGhost ){
950	>
951	>	exI = theBends[j].a;
952	>	exJ = theBends[j].c;
953	>	}
954	>	else{
955
956	<	#ifdef IS_MPI
956	>	exI = theBends[j].a;
957	>	exJ = theBends[j].b;
958		}
937	–	#endif //is_mpi
959
960	<	molIndex++;
961	<	}
962	<	}
960	>	// exclude_I must always be the smaller of the pair
961	>	if( exI > exJ ){
962	>	tempEx = exI;
963	>	exI = exJ;
964	>	exJ = tempEx;
965	>	}
966	>	#ifdef IS_MPI
967	>	tempEx = exI;
968	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
969	>	tempEx = exJ;
970	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
971	>
972	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
973	>	#else // isn't MPI
974	>	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
975	>	#endif //is_mpi
976	>	}
977	>	excludeOffset += info.nBends;
978
979	<	the_ff->initializeBonds( the_bonds );
980	<	}
979	>	for(j=0; j<info.nTorsions; j++){
980	>
981	>	currentTorsion = comp_stamps[stampID]->getTorsion( j );
982	>	theTorsions[j].a = currentTorsion->getA() + atomOffset;
983	>	theTorsions[j].b = currentTorsion->getB() + atomOffset;
984	>	theTorsions[j].c = currentTorsion->getC() + atomOffset;
985	>	theTorsions[j].d = currentTorsion->getD() + atomOffset;
986	>
987	>	exI = theTorsions[j].a;
988	>	exJ = theTorsions[j].d;
989
990	<	void SimSetup::makeBends( void ){
990	>	// exclude_I must always be the smaller of the pair
991	>	if( exI > exJ ){
992	>	tempEx = exI;
993	>	exI = exJ;
994	>	exJ = tempEx;
995	>	}
996	>	#ifdef IS_MPI
997	>	tempEx = exI;
998	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
999	>	tempEx = exJ;
1000	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
1001	>
1002	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
1003	>	#else // isn't MPI
1004	>	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
1005	>	#endif //is_mpi
1006	>	}
1007	>	excludeOffset += info.nTorsions;
1008
1009	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1010	<	bend_set* the_bends;
950	<	BendStamp* current_bend;
951	<	LinkedAssign* extras;
952	<	LinkedAssign* current_extra;
953	<
1009	>
1010	>	// send the arrays off to the forceField for init.
1011
1012	<	the_bends = new bend_set[tot_bends];
1013	<	index = 0;
1014	<	offset = 0;
1015	<	molIndex = 0;
959	<	for( i=0; i<n_components; i++ ){
1012	>	the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
1013	>	the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
1014	>	the_ff->initializeBends( info.nBends, info.myBends, theBends );
1015	>	the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
1016
961	–	for( j=0; j<components_nmol[i]; j++ ){
1017
1018	<	#ifdef IS_MPI
964	<	if( mpiSim->getMyMolStart() <= molIndex &&
965	<	molIndex <= mpiSim->getMyMolEnd() ){
966	<	#endif // is_mpi
967	<
968	<	for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
969	<
970	<	current_bend = comp_stamps[i]->getBend( k );
971	<	the_bends[index].a = current_bend->getA() + offset;
972	<	the_bends[index].b = current_bend->getB() + offset;
973	<	the_bends[index].c = current_bend->getC() + offset;
974	<
975	<	if( current_bend->haveExtras() ){
976	<
977	<	extras = current_bend->getExtras();
978	<	current_extra = extras;
979	<
980	<	while( current_extra != NULL ){
981	<	if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
982	<
983	<	switch( current_extra->getType() ){
984	<
985	<	case 0:
986	<	the_bends[index].ghost =
987	<	current_extra->getInt() + offset;
988	<	the_bends[index].isGhost = 1;
989	<	break;
990	<
991	<	case 1:
992	<	the_bends[index].ghost =
993	<	(int)current_extra->getDouble() + offset;
994	<	the_bends[index].isGhost = 1;
995	<	break;
996	<
997	<	default:
998	<	sprintf( painCave.errMsg,
999	<	"SimSetup Error: ghostVectorSource was neiter a "
1000	<	"double nor an int.\n"
1001	<	"-->Bend[%d] in %s\n",
1002	<	k, comp_stamps[i]->getID() );
1003	<	painCave.isFatal = 1;
1004	<	simError();
1005	<	}
1006	<	}
1007	<
1008	<	else{
1009	<
1010	<	sprintf( painCave.errMsg,
1011	<	"SimSetup Error: unhandled bend assignment:\n"
1012	<	" -->%s in Bend[%d] in %s\n",
1013	<	current_extra->getlhs(),
1014	<	k, comp_stamps[i]->getID() );
1015	<	painCave.isFatal = 1;
1016	<	simError();
1017	<	}
1018	<
1019	<	current_extra = current_extra->getNext();
1020	<	}
1021	<	}
1022	<
1023	<	if( !the_bends[index].isGhost ){
1024	<
1025	<	exI = the_bends[index].a;
1026	<	exJ = the_bends[index].c;
1027	<	}
1028	<	else{
1029	<
1030	<	exI = the_bends[index].a;
1031	<	exJ = the_bends[index].b;
1032	<	}
1033	<
1034	<	// exclude_I must always be the smaller of the pair
1035	<	if( exI > exJ ){
1036	<	tempEx = exI;
1037	<	exI = exJ;
1038	<	exJ = tempEx;
1039	<	}
1018	>	the_molecules[i].initialize( info );
1019
1041	–
1042	–	#ifdef IS_MPI
1043	–
1044	–	the_excludes[(index + tot_bonds)*2] =
1045	–	the_atoms[exI]->getGlobalIndex() + 1;
1046	–	the_excludes[(index + tot_bonds)*2 + 1] =
1047	–	the_atoms[exJ]->getGlobalIndex() + 1;
1048	–
1049	–	#else // isn't MPI
1050	–
1051	–	the_excludes[(index + tot_bonds)*2] = exI + 1;
1052	–	the_excludes[(index + tot_bonds)*2 + 1] = exJ + 1;
1053	–	// fortran index from 1 (hence the +1 in the indexing)
1054	–	#endif //is_mpi
1055	–
1056	–
1057	–	// increment the index and repeat;
1058	–	index++;
1059	–	}
1060	–	offset += comp_stamps[i]->getNAtoms();
1061	–
1062	–	#ifdef IS_MPI
1063	–	}
1064	–	#endif //is_mpi
1020
1021	<	molIndex++;
1022	<	}
1021	>	atomOffset += info.nAtoms;
1022	>	delete[] theBonds;
1023	>	delete[] theBends;
1024	>	delete[] theTorsions;
1025		}
1026
1027		#ifdef IS_MPI
1028	<	sprintf( checkPointMsg,
1072	<	"Successfully created the bends list.\n" );
1028	>	sprintf( checkPointMsg, "all molecules initialized succesfully" );
1029		MPIcheckPoint();
1030		#endif // is_mpi
1075	–
1031
1032	<	the_ff->initializeBends( the_bends );
1033	<	}
1032	>	// clean up the forcefield
1033	>	the_ff->calcRcut();
1034	>	the_ff->cleanMe();
1035
1080	–	void SimSetup::makeTorsions( void ){
1081	–
1082	–	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1083	–	torsion_set* the_torsions;
1084	–	TorsionStamp* current_torsion;
1085	–
1086	–	the_torsions = new torsion_set[tot_torsions];
1087	–	index = 0;
1088	–	offset = 0;
1089	–	molIndex = 0;
1090	–	for( i=0; i<n_components; i++ ){
1091	–
1092	–	for( j=0; j<components_nmol[i]; j++ ){
1093	–
1094	–	#ifdef IS_MPI
1095	–	if( mpiSim->getMyMolStart() <= molIndex &&
1096	–	molIndex <= mpiSim->getMyMolEnd() ){
1097	–	#endif // is_mpi
1098	–
1099	–	for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
1100	–
1101	–	current_torsion = comp_stamps[i]->getTorsion( k );
1102	–	the_torsions[index].a = current_torsion->getA() + offset;
1103	–	the_torsions[index].b = current_torsion->getB() + offset;
1104	–	the_torsions[index].c = current_torsion->getC() + offset;
1105	–	the_torsions[index].d = current_torsion->getD() + offset;
1106	–
1107	–	exI = the_torsions[index].a;
1108	–	exJ = the_torsions[index].d;
1109	–
1110	–
1111	–	// exclude_I must always be the smaller of the pair
1112	–	if( exI > exJ ){
1113	–	tempEx = exI;
1114	–	exI = exJ;
1115	–	exJ = tempEx;
1116	–	}
1117	–
1118	–
1119	–	#ifdef IS_MPI
1120	–
1121	–	the_excludes[(index + tot_bonds + tot_bends)*2] =
1122	–	the_atoms[exI]->getGlobalIndex() + 1;
1123	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] =
1124	–	the_atoms[exJ]->getGlobalIndex() + 1;
1125	–
1126	–	#else // isn't MPI
1127	–
1128	–	the_excludes[(index + tot_bonds + tot_bends)*2] = exI + 1;
1129	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] = exJ + 1;
1130	–	// fortran indexes from 1 (hence the +1 in the indexing)
1131	–	#endif //is_mpi
1132	–
1133	–
1134	–	// increment the index and repeat;
1135	–	index++;
1136	–	}
1137	–	offset += comp_stamps[i]->getNAtoms();
1138	–
1139	–	#ifdef IS_MPI
1140	–	}
1141	–	#endif //is_mpi
1142	–
1143	–	molIndex++;
1144	–	}
1145	–	}
1146	–
1147	–	the_ff->initializeTorsions( the_torsions );
1036		}
1037
1038		void SimSetup::initFromBass( void ){

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents): Revision 407 by mmeineke, Wed Mar 26 20:22:02 2003 UTC vs. Revision 489 by mmeineke, Thu Apr 10 20:08:56 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 407 by mmeineke, Wed Mar 26 20:22:02 2003 UTC vs.
Revision 489 by mmeineke, Thu Apr 10 20:08:56 2003 UTC