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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 1032 by tim, Fri Feb 6 19:05:07 2004 UTC vs.
Revision 1129 by tim, Thu Apr 22 03:29:30 2004 UTC

#	Line 9 \| Line 9
9		#include "parse_me.h"
10		#include "Integrator.hpp"
11		#include "simError.h"
12	<	#include "ConjugateMinimizer.hpp"
12	>	#include "RigidBody.hpp"
13	>	//#include "ConjugateMinimizer.hpp"
14	>	#include "OOPSEMinimizer.hpp"
15
16		#ifdef IS_MPI
17		#include "mpiBASS.h"
#	Line 28 \| Line 30
30		#define FF_DUFF 0
31		#define FF_LJ 1
32		#define FF_EAM 2
33	<	#define FF_H2O 3
33	>	#define FF_H2O 3
34
35		using namespace std;
36
#	Line 145 \| Line 147 \| void SimSetup::createSim(void){
147		// make the output filenames
148
149		makeOutNames();
148	–
149	–	if (globals->haveMinimizer())
150	–	// make minimizer
151	–	makeMinimizer();
152	–	else
153	–	// make the integrator
154	–	makeIntegrator();
150
151		#ifdef IS_MPI
152		mpiSim->mpiRefresh();
#	Line 160 \| Line 155 \| void SimSetup::createSim(void){
155		// initialize the Fortran
156
157		initFortran();
158	+
159	+	if (globals->haveMinimizer())
160	+	// make minimizer
161	+	makeMinimizer();
162	+	else
163	+	// make the integrator
164	+	makeIntegrator();
165	+
166		}
167
168
169		void SimSetup::makeMolecules(void){
170	<	int k;
171	<	int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
170	>	int i, j, k;
171	>	int exI, exJ, exK, exL, slI, slJ;
172	>	int tempI, tempJ, tempK, tempL;
173	>	int molI;
174	>	int stampID, atomOffset, rbOffset;
175		molInit molInfo;
176		DirectionalAtom* dAtom;
177	+	RigidBody* myRB;
178	+	StuntDouble* mySD;
179		LinkedAssign* extras;
180		LinkedAssign* current_extra;
181		AtomStamp* currentAtom;
182		BondStamp* currentBond;
183		BendStamp* currentBend;
184		TorsionStamp* currentTorsion;
185	+	RigidBodyStamp* currentRigidBody;
186
187		bond_pair* theBonds;
188		bend_set* theBends;
189		torsion_set* theTorsions;
190
191	+	set<int> skipList;
192	+
193	+	double phi, theta, psi;
194	+	char* molName;
195	+	char rbName[100];
196	+
197		//init the forceField paramters
198
199		the_ff->readParams();
200
186	–
201		// init the atoms
202
203	<	double phi, theta, psi;
190	<	double sux, suy, suz;
191	<	double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
192	<	double ux, uy, uz, u, uSqr;
203	>	int nMembers, nNew, rb1, rb2;
204
205		for (k = 0; k < nInfo; k++){
206		the_ff->setSimInfo(&(info[k]));
207
208		atomOffset = 0;
209	<	excludeOffset = 0;
209	>
210		for (i = 0; i < info[k].n_mol; i++){
211		stampID = info[k].molecules[i].getStampID();
212	+	molName = comp_stamps[stampID]->getID();
213
214		molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
215		molInfo.nBonds = comp_stamps[stampID]->getNBonds();
216		molInfo.nBends = comp_stamps[stampID]->getNBends();
217		molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
218	<	molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
219	<
218	>	molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
219	>
220		molInfo.myAtoms = &(info[k].atoms[atomOffset]);
209	–	molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
210	–	molInfo.myBonds = new Bond * [molInfo.nBonds];
211	–	molInfo.myBends = new Bend * [molInfo.nBends];
212	–	molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
221
222	+	if (molInfo.nBonds > 0)
223	+	molInfo.myBonds = new (Bond *) [molInfo.nBonds];
224	+	else
225	+	molInfo.myBonds = NULL;
226	+
227	+	if (molInfo.nBends > 0)
228	+	molInfo.myBends = new (Bend *) [molInfo.nBends];
229	+	else
230	+	molInfo.myBends = NULL;
231	+
232	+	if (molInfo.nTorsions > 0)
233	+	molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
234	+	else
235	+	molInfo.myTorsions = NULL;
236	+
237		theBonds = new bond_pair[molInfo.nBonds];
238		theBends = new bend_set[molInfo.nBends];
239		theTorsions = new torsion_set[molInfo.nTorsions];
240	<
240	>
241		// make the Atoms
242
243		for (j = 0; j < molInfo.nAtoms; j++){
244		currentAtom = comp_stamps[stampID]->getAtom(j);
245	+
246		if (currentAtom->haveOrientation()){
247		dAtom = new DirectionalAtom((j + atomOffset),
248		info[k].getConfiguration());
#	Line 232 \| Line 256 \| void SimSetup::makeMolecules(void){
256		phi = currentAtom->getEulerPhi() * M_PI / 180.0;
257		theta = currentAtom->getEulerTheta() * M_PI / 180.0;
258		psi = currentAtom->getEulerPsi()* M_PI / 180.0;
259	+
260	+	dAtom->setUnitFrameFromEuler(phi, theta, psi);
261
262	<	Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
263	<	Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
238	<	Axz = sin(theta) * sin(psi);
239	<
240	<	Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
241	<	Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
242	<	Ayz = sin(theta) * cos(psi);
243	<
244	<	Azx = sin(phi) * sin(theta);
245	<	Azy = -cos(phi) * sin(theta);
246	<	Azz = cos(theta);
262	>	}
263	>	else{
264
265	<	sux = 0.0;
249	<	suy = 0.0;
250	<	suz = 1.0;
265	>	molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
266
252	–	ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
253	–	uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
254	–	uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
255	–
256	–	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
257	–
258	–	u = sqrt(uSqr);
259	–	ux = ux / u;
260	–	uy = uy / u;
261	–	uz = uz / u;
262	–
263	–	dAtom->setSUx(ux);
264	–	dAtom->setSUy(uy);
265	–	dAtom->setSUz(uz);
267		}
267	–	else{
268	–	molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
269	–	info[k].getConfiguration());
270	–	}
271	–	molInfo.myAtoms[j]->setType(currentAtom->getType());
268
269	+	molInfo.myAtoms[j]->setType(currentAtom->getType());
270		#ifdef IS_MPI
271
272	<	molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
272	>	molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
273
274		#endif // is_mpi
275		}
#	Line 283 \| Line 280 \| void SimSetup::makeMolecules(void){
280		theBonds[j].a = currentBond->getA() + atomOffset;
281		theBonds[j].b = currentBond->getB() + atomOffset;
282
283	<	exI = theBonds[j].a;
284	<	exJ = theBonds[j].b;
283	>	tempI = theBonds[j].a;
284	>	tempJ = theBonds[j].b;
285
289	–	// exclude_I must always be the smaller of the pair
290	–	if (exI > exJ){
291	–	tempEx = exI;
292	–	exI = exJ;
293	–	exJ = tempEx;
294	–	}
286		#ifdef IS_MPI
287	<	tempEx = exI;
288	<	exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
289	<	tempEx = exJ;
290	<	exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
291	<
292	<	info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
302	<	#else // isn't MPI
287	>	exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
288	>	exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
289	>	#else
290	>	exI = tempI + 1;
291	>	exJ = tempJ + 1;
292	>	#endif
293
294	<	info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
305	<	#endif //is_mpi
294	>	info[k].excludes->addPair(exI, exJ);
295		}
307	–	excludeOffset += molInfo.nBonds;
296
297		//make the bends
298		for (j = 0; j < molInfo.nBends; j++){
#	Line 354 \| Line 342 \| void SimSetup::makeMolecules(void){
342		}
343		}
344
345	<	if (!theBends[j].isGhost){
346	<	exI = theBends[j].a;
347	<	exJ = theBends[j].c;
348	<	}
349	<	else{
362	<	exI = theBends[j].a;
363	<	exJ = theBends[j].b;
364	<	}
365	<
366	<	// exclude_I must always be the smaller of the pair
367	<	if (exI > exJ){
368	<	tempEx = exI;
369	<	exI = exJ;
370	<	exJ = tempEx;
371	<	}
345	>	if (theBends[j].isGhost) {
346	>
347	>	tempI = theBends[j].a;
348	>	tempJ = theBends[j].b;
349	>
350		#ifdef IS_MPI
351	<	tempEx = exI;
352	<	exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
353	<	tempEx = exJ;
354	<	exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
351	>	exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
352	>	exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
353	>	#else
354	>	exI = tempI + 1;
355	>	exJ = tempJ + 1;
356	>	#endif
357	>	info[k].excludes->addPair(exI, exJ);
358
359	<	info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
360	<	#else // isn't MPI
361	<	info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
362	<	#endif //is_mpi
359	>	} else {
360	>
361	>	tempI = theBends[j].a;
362	>	tempJ = theBends[j].b;
363	>	tempK = theBends[j].c;
364	>
365	>	#ifdef IS_MPI
366	>	exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
367	>	exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
368	>	exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
369	>	#else
370	>	exI = tempI + 1;
371	>	exJ = tempJ + 1;
372	>	exK = tempK + 1;
373	>	#endif
374	>
375	>	info[k].excludes->addPair(exI, exK);
376	>	info[k].excludes->addPair(exI, exJ);
377	>	info[k].excludes->addPair(exJ, exK);
378	>	}
379		}
383	–	excludeOffset += molInfo.nBends;
380
381		for (j = 0; j < molInfo.nTorsions; j++){
382		currentTorsion = comp_stamps[stampID]->getTorsion(j);
#	Line 389 \| Line 385 \| void SimSetup::makeMolecules(void){
385		theTorsions[j].c = currentTorsion->getC() + atomOffset;
386		theTorsions[j].d = currentTorsion->getD() + atomOffset;
387
388	<	exI = theTorsions[j].a;
389	<	exJ = theTorsions[j].d;
388	>	tempI = theTorsions[j].a;
389	>	tempJ = theTorsions[j].b;
390	>	tempK = theTorsions[j].c;
391	>	tempL = theTorsions[j].d;
392
395	–	// exclude_I must always be the smaller of the pair
396	–	if (exI > exJ){
397	–	tempEx = exI;
398	–	exI = exJ;
399	–	exJ = tempEx;
400	–	}
393		#ifdef IS_MPI
394	<	tempEx = exI;
395	<	exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
396	<	tempEx = exJ;
397	<	exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
394	>	exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
395	>	exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
396	>	exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
397	>	exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
398	>	#else
399	>	exI = tempI + 1;
400	>	exJ = tempJ + 1;
401	>	exK = tempK + 1;
402	>	exL = tempL + 1;
403	>	#endif
404
405	<	info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
406	<	#else // isn't MPI
407	<	info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
408	<	#endif //is_mpi
405	>	info[k].excludes->addPair(exI, exJ);
406	>	info[k].excludes->addPair(exI, exK);
407	>	info[k].excludes->addPair(exI, exL);
408	>	info[k].excludes->addPair(exJ, exK);
409	>	info[k].excludes->addPair(exJ, exL);
410	>	info[k].excludes->addPair(exK, exL);
411		}
412	–	excludeOffset += molInfo.nTorsions;
412
413	+
414	+	molInfo.myRigidBodies.clear();
415	+
416	+	for (j = 0; j < molInfo.nRigidBodies; j++){
417
418	<	// send the arrays off to the forceField for init.
418	>	currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
419	>	nMembers = currentRigidBody->getNMembers();
420	>
421	>	// Create the Rigid Body:
422	>
423	>	myRB = new RigidBody();
424	>
425	>	sprintf(rbName,"%s_RB_%d", molName, j);
426	>	myRB->setType(rbName);
427	>
428	>	for (rb1 = 0; rb1 < nMembers; rb1++) {
429	>
430	>	// molI is atom numbering inside this molecule
431	>	molI = currentRigidBody->getMember(rb1);
432	>
433	>	// tempI is atom numbering on local processor
434	>	tempI = molI + atomOffset;
435	>
436	>	// currentAtom is the AtomStamp (which we need for
437	>	// rigid body reference positions)
438	>	currentAtom = comp_stamps[stampID]->getAtom(molI);
439	>
440	>	// When we add to the rigid body, add the atom itself and
441	>	// the stamp info:
442	>
443	>	myRB->addAtom(info[k].atoms[tempI], currentAtom);
444	>
445	>	// Add this atom to the Skip List for the integrators
446	>	#ifdef IS_MPI
447	>	slI = info[k].atoms[tempI]->getGlobalIndex();
448	>	#else
449	>	slI = tempI;
450	>	#endif
451	>	skipList.insert(slI);
452	>
453	>	}
454	>
455	>	for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
456	>	for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
457	>
458	>	tempI = currentRigidBody->getMember(rb1);
459	>	tempJ = currentRigidBody->getMember(rb2);
460	>
461	>	// Some explanation is required here.
462	>	// Fortran indexing starts at 1, while c indexing starts at 0
463	>	// Also, in parallel computations, the GlobalIndex is
464	>	// used for the exclude list:
465	>
466	>	#ifdef IS_MPI
467	>	exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
468	>	exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
469	>	#else
470	>	exI = molInfo.myAtoms[tempI]->getIndex() + 1;
471	>	exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
472	>	#endif
473	>
474	>	info[k].excludes->addPair(exI, exJ);
475	>
476	>	}
477	>	}
478	>
479	>	molInfo.myRigidBodies.push_back(myRB);
480	>	info[k].rigidBodies.push_back(myRB);
481	>	}
482	>
483	>
484	>	// After this is all set up, scan through the atoms to
485	>	// see if they can be added to the integrableObjects:
486	>
487	>	molInfo.myIntegrableObjects.clear();
488	>
489	>
490	>	for (j = 0; j < molInfo.nAtoms; j++){
491	>
492	>	#ifdef IS_MPI
493	>	slJ = molInfo.myAtoms[j]->getGlobalIndex();
494	>	#else
495	>	slJ = j+atomOffset;
496	>	#endif
497	>
498	>	// if they aren't on the skip list, then they can be integrated
499
500	+	if (skipList.find(slJ) == skipList.end()) {
501	+	mySD = (StuntDouble *) molInfo.myAtoms[j];
502	+	info[k].integrableObjects.push_back(mySD);
503	+	molInfo.myIntegrableObjects.push_back(mySD);
504	+	}
505	+	}
506	+
507	+	// all rigid bodies are integrated:
508	+
509	+	for (j = 0; j < molInfo.nRigidBodies; j++) {
510	+	mySD = (StuntDouble *) molInfo.myRigidBodies[j];
511	+	info[k].integrableObjects.push_back(mySD);
512	+	molInfo.myIntegrableObjects.push_back(mySD);
513	+	}
514	+
515	+
516	+	// send the arrays off to the forceField for init.
517	+
518		the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
519		the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
520		the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
521		the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
522		theTorsions);
523
423	–
524		info[k].molecules[i].initialize(molInfo);
525
526
#	Line 428 \| Line 528 \| void SimSetup::makeMolecules(void){
528		delete[] theBonds;
529		delete[] theBends;
530		delete[] theTorsions;
531	<	}
531	>	}
532		}
533
534		#ifdef IS_MPI
#	Line 438 \| Line 538 \| void SimSetup::makeMolecules(void){
538
539		// clean up the forcefield
540
541	<	the_ff->calcRcut();
541	>	if (!globals->haveLJrcut()){
542	>
543	>	the_ff->calcRcut();
544	>
545	>	} else {
546	>
547	>	the_ff->setRcut( globals->getLJrcut() );
548	>	}
549	>
550		the_ff->cleanMe();
551		}
552
#	Line 726 \| Line 834 \| void SimSetup::gatherInfo(void){
834		}
835
836		//check whether sample time, status time, thermal time and reset time are divisble by dt
837	<	if (!isDivisible(globals->getSampleTime(), globals->getDt())){
837	>	if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
838		sprintf(painCave.errMsg,
839		"Sample time is not divisible by dt.\n"
840		"\tThis will result in samples that are not uniformly\n"
#	Line 736 \| Line 844 \| void SimSetup::gatherInfo(void){
844		simError();
845		}
846
847	<	if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
847	>	if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
848		sprintf(painCave.errMsg,
849		"Status time is not divisible by dt.\n"
850		"\tThis will result in status reports that are not uniformly\n"
#	Line 772 \| Line 880 \| void SimSetup::gatherInfo(void){
880		if (globals->haveSampleTime()){
881		info[i].sampleTime = globals->getSampleTime();
882		info[i].statusTime = info[i].sampleTime;
775	–	info[i].thermalTime = info[i].sampleTime;
883		}
884		else{
885		info[i].sampleTime = globals->getRunTime();
886		info[i].statusTime = info[i].sampleTime;
780	–	info[i].thermalTime = info[i].sampleTime;
887		}
888
889		if (globals->haveStatusTime()){
#	Line 786 \| Line 892 \| void SimSetup::gatherInfo(void){
892
893		if (globals->haveThermalTime()){
894		info[i].thermalTime = globals->getThermalTime();
895	+	} else {
896	+	info[i].thermalTime = globals->getRunTime();
897		}
898
899		info[i].resetIntegrator = 0;
#	Line 856 \| Line 964 \| void SimSetup::finalInfoCheck(void){
964		void SimSetup::finalInfoCheck(void){
965		int index;
966		int usesDipoles;
967	+	int usesCharges;
968		int i;
969
970		for (i = 0; i < nInfo; i++){
#	Line 867 \| Line 976 \| void SimSetup::finalInfoCheck(void){
976		usesDipoles = (info[i].atoms[index])->hasDipole();
977		index++;
978		}
979	<
979	>	index = 0;
980	>	usesCharges = 0;
981	>	while ((index < info[i].n_atoms) && !usesCharges){
982	>	usesCharges= (info[i].atoms[index])->hasCharge();
983	>	index++;
984	>	}
985		#ifdef IS_MPI
986		int myUse = usesDipoles;
987		MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
#	Line 918 \| Line 1032 \| void SimSetup::finalInfoCheck(void){
1032		info[i].dielectric = globals->getDielectric();
1033		}
1034		else{
1035	<	if (usesDipoles){
1035	>	if (usesDipoles \|\| usesCharges){
1036		if (!globals->haveECR()){
1037		sprintf(painCave.errMsg,
1038		"SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
#	Line 1220 \| Line 1334 \| void SimSetup::calcSysValues(void){
1334		tot_bonds = 0;
1335		tot_bends = 0;
1336		tot_torsions = 0;
1337	+	tot_rigid = 0;
1338		for (i = 0; i < n_components; i++){
1339		tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1340		tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1341		tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1342		tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1343	+	tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1344		}
1345	<
1345	>
1346		tot_SRI = tot_bonds + tot_bends + tot_torsions;
1347		molMembershipArray = new int[tot_atoms];
1348
#	Line 1248 \| Line 1364 \| void SimSetup::mpiMolDivide(void){
1364		int i, j, k;
1365		int localMol, allMol;
1366		int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1367	+	int local_rigid;
1368	+	vector<int> globalMolIndex;
1369
1370		mpiSim = new mpiSimulation(info);
1371
1372	<	globalIndex = mpiSim->divideLabor();
1372	>	mpiSim->divideLabor();
1373	>	globalAtomIndex = mpiSim->getGlobalAtomIndex();
1374	>	//globalMolIndex = mpiSim->getGlobalMolIndex();
1375
1376		// set up the local variables
1377
#	Line 1264 \| Line 1384 \| void SimSetup::mpiMolDivide(void){
1384		local_bonds = 0;
1385		local_bends = 0;
1386		local_torsions = 0;
1387	<	globalAtomIndex = 0;
1388	<
1387	>	local_rigid = 0;
1388	>	globalAtomCounter = 0;
1389
1390		for (i = 0; i < n_components; i++){
1391		for (j = 0; j < components_nmol[i]; j++){
#	Line 1274 \| Line 1394 \| void SimSetup::mpiMolDivide(void){
1394		local_bonds += comp_stamps[i]->getNBonds();
1395		local_bends += comp_stamps[i]->getNBends();
1396		local_torsions += comp_stamps[i]->getNTorsions();
1397	+	local_rigid += comp_stamps[i]->getNRigidBodies();
1398		localMol++;
1399		}
1400		for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1401	<	info[0].molMembershipArray[globalAtomIndex] = allMol;
1402	<	globalAtomIndex++;
1401	>	info[0].molMembershipArray[globalAtomCounter] = allMol;
1402	>	globalAtomCounter++;
1403		}
1404
1405		allMol++;
#	Line 1287 \| Line 1408 \| void SimSetup::mpiMolDivide(void){
1408		local_SRI = local_bonds + local_bends + local_torsions;
1409
1410		info[0].n_atoms = mpiSim->getMyNlocal();
1411	+
1412
1413		if (local_atoms != info[0].n_atoms){
1414		sprintf(painCave.errMsg,
#	Line 1319 \| Line 1441 \| void SimSetup::makeSysArrays(void){
1441
1442		Atom** the_atoms;
1443		Molecule* the_molecules;
1322	–	Exclude** the_excludes;
1444
1324	–
1445		for (l = 0; l < nInfo; l++){
1446		// create the atom and short range interaction arrays
1447
#	Line 1347 \| Line 1467 \| void SimSetup::makeSysArrays(void){
1467		#else // is_mpi
1468
1469		molIndex = 0;
1470	<	globalAtomIndex = 0;
1470	>	globalAtomCounter = 0;
1471		for (i = 0; i < n_components; i++){
1472		for (j = 0; j < components_nmol[i]; j++){
1473		the_molecules[molIndex].setStampID(i);
1474		the_molecules[molIndex].setMyIndex(molIndex);
1475		the_molecules[molIndex].setGlobalIndex(molIndex);
1476		for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1477	<	info[l].molMembershipArray[globalAtomIndex] = molIndex;
1478	<	globalAtomIndex++;
1477	>	info[l].molMembershipArray[globalAtomCounter] = molIndex;
1478	>	globalAtomCounter++;
1479		}
1480		molIndex++;
1481		}
#	Line 1364 \| Line 1484 \| void SimSetup::makeSysArrays(void){
1484
1485		#endif // is_mpi
1486
1487	<
1488	<	if (info[l].n_SRI){
1489	<	Exclude::createArray(info[l].n_SRI);
1370	<	the_excludes = new Exclude * [info[l].n_SRI];
1371	<	for (int ex = 0; ex < info[l].n_SRI; ex++){
1372	<	the_excludes[ex] = new Exclude(ex);
1373	<	}
1374	<	info[l].globalExcludes = new int;
1375	<	info[l].n_exclude = info[l].n_SRI;
1376	<	}
1377	<	else{
1378	<	Exclude::createArray(1);
1379	<	the_excludes = new Exclude * ;
1380	<	the_excludes[0] = new Exclude(0);
1381	<	the_excludes[0]->setPair(0, 0);
1382	<	info[l].globalExcludes = new int;
1383	<	info[l].globalExcludes[0] = 0;
1384	<	info[l].n_exclude = 0;
1385	<	}
1386	<
1487	>	info[l].globalExcludes = new int;
1488	>	info[l].globalExcludes[0] = 0;
1489	>
1490		// set the arrays into the SimInfo object
1491
1492		info[l].atoms = the_atoms;
1493		info[l].molecules = the_molecules;
1494		info[l].nGlobalExcludes = 0;
1392	–	info[l].excludes = the_excludes;
1495
1496		the_ff->setSimInfo(info);
1497		}
#	Line 1662 \| Line 1764 \| void SimSetup::setupZConstraint(SimInfo& theInfo){
1764		}
1765
1766		theInfo.addProperty(zconsForcePolicy);
1767	+
1768	+	//set zcons gap
1769	+	DoubleData* zconsGap = new DoubleData();
1770	+	zconsGap->setID(ZCONSGAP_ID);
1771	+
1772	+	if (globals->haveZConsGap()){
1773	+	zconsGap->setData(globals->getZconsGap());
1774	+	theInfo.addProperty(zconsGap);
1775	+	}
1776
1777	+	//set zcons fixtime
1778	+	DoubleData* zconsFixtime = new DoubleData();
1779	+	zconsFixtime->setID(ZCONSFIXTIME_ID);
1780	+
1781	+	if (globals->haveZConsFixTime()){
1782	+	zconsFixtime->setData(globals->getZconsFixtime());
1783	+	theInfo.addProperty(zconsFixtime);
1784	+	}
1785	+
1786	+	//set zconsUsingSMD
1787	+	IntData* zconsUsingSMD = new IntData();
1788	+	zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1789	+
1790	+	if (globals->haveZConsUsingSMD()){
1791	+	zconsUsingSMD->setData(globals->getZconsUsingSMD());
1792	+	theInfo.addProperty(zconsUsingSMD);
1793	+	}
1794	+
1795		//Determine the name of ouput file and add it into SimInfo's property list
1796		//Be careful, do not use inFileName, since it is a pointer which
1797		//point to a string at master node, and slave nodes do not contain that string
#	Line 1692 \| Line 1821 \| void SimSetup::setupZConstraint(SimInfo& theInfo){
1821		tempParaItem.zPos = zconStamp[i]->getZpos();
1822		tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1823		tempParaItem.kRatio = zconStamp[i]->getKratio();
1824	<
1824	>	tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1825	>	tempParaItem.cantVel = zconStamp[i]->getCantVel();
1826		zconsParaData->addItem(tempParaItem);
1827		}
1828
#	Line 1713 \| Line 1843 \| void SimSetup::makeMinimizer(){
1843
1844		void SimSetup::makeMinimizer(){
1845
1846	<	OOPSEMinimizerBase* myOOPSEMinimizerBase;
1717	<	ObjFunctor1 * objFunc;
1718	<	OutputFunctor* outputFunc;
1719	<	ConcreteNLModel1* nlp;
1846	>	OOPSEMinimizer* myOOPSEMinimizer;
1847		MinimizerParameterSet* param;
1848	<	ConjugateMinimizerBase* minimizer;
1722	<	int dim;
1848	>	char minimizerName[100];
1849
1850		for (int i = 0; i < nInfo; i++){
1851	<	//creat
1726	<	myOOPSEMinimizerBase = new OOPSEMinimizerBase(&(info[i]), the_ff);
1727	<
1728	<	//creat the object functor;
1729	<	objFunc = (ObjFunctor1*) new ClassMemObjFunctor1<OOPSEMinimizerBase>
1730	<	(myOOPSEMinimizerBase, &OOPSEMinimizerBase::calcGradient);
1731	<
1732	<	//creat output functor;
1733	<	outputFunc = new ClassMemOutputFunctor<OOPSEMinimizerBase>
1734	<	(myOOPSEMinimizerBase, &OOPSEMinimizerBase::output);
1735	<
1736	<	//creat nonlinear model
1737	<	dim = myOOPSEMinimizerBase->getDim();
1738	<	nlp = new ConcreteNLModel1(dim, objFunc);
1739	<
1851	>
1852		//prepare parameter set for minimizer
1853		param = new MinimizerParameterSet();
1854		param->setDefaultParameter();
#	Line 1761 \| Line 1873 \| void SimSetup::makeMinimizer(){
1873		param->setWriteFrq(globals->getMinWriteFrq());
1874		}
1875
1876	<	if (globals->haveMinResetFrq()){
1877	<	param->setResetFrq(globals->getMinResetFrq());
1876	>	if (globals->haveMinStepSize()){
1877	>	param->setStepSize(globals->getMinStepSize());
1878		}
1879
1880		if (globals->haveMinLSMaxIter()){
#	Line 1772 \| Line 1884 \| void SimSetup::makeMinimizer(){
1884		if (globals->haveMinLSTol()){
1885		param->setLineSearchTol(globals->getMinLSTol());
1886		}
1775	–
1776	–	//creat the minimizer
1777	–	minimizer = new PRCGMinimizer(nlp, param);
1778	–	minimizer->setLineSearchStrategy(nlp, GoldenSection);
1779	–	minimizer->setOutputFunctor(outputFunc);
1887
1888	+	strcpy(minimizerName, globals->getMinimizer());
1889	+
1890	+	if (!strcasecmp(minimizerName, "CG")){
1891	+	myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1892	+	}
1893	+	else if (!strcasecmp(minimizerName, "SD")){
1894	+	//myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
1895	+	myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
1896	+	}
1897	+	else{
1898	+	sprintf(painCave.errMsg,
1899	+	"SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
1900	+	painCave.isFatal = 0;
1901	+	simError();
1902	+
1903	+	myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
1904	+	}
1905	+	info[i].the_integrator = myOOPSEMinimizer;
1906	+
1907		//store the minimizer into simInfo
1908	<	info[i].the_minimizer = minimizer;
1908	>	info[i].the_minimizer = myOOPSEMinimizer;
1909		info[i].has_minimizer = true;
1910		}
1911

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents): Revision 1032 by tim, Fri Feb 6 19:05:07 2004 UTC vs. Revision 1129 by tim, Thu Apr 22 03:29:30 2004 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 1032 by tim, Fri Feb 6 19:05:07 2004 UTC vs.
Revision 1129 by tim, Thu Apr 22 03:29:30 2004 UTC