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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 1096 by tim, Wed Mar 17 14:22:59 2004 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

#	Line 9 \| Line 9
9		#include "parse_me.h"
10		#include "Integrator.hpp"
11		#include "simError.h"
12	+	#include "RigidBody.hpp"
13		//#include "ConjugateMinimizer.hpp"
14		#include "OOPSEMinimizer.hpp"
15
#	Line 165 \| Line 166 \| void SimSetup::makeMolecules(void){
166
167
168		void SimSetup::makeMolecules(void){
169	<	int k;
170	<	int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
169	>	int i, j, k;
170	>	int exI, exJ, exK, exL, slI;
171	>	int tempI, tempJ, tempK, tempL;
172	>	int molI;
173	>	int stampID, atomOffset, rbOffset;
174		molInit molInfo;
175		DirectionalAtom* dAtom;
176	+	RigidBody* myRB;
177	+	StuntDouble* mySD;
178		LinkedAssign* extras;
179		LinkedAssign* current_extra;
180		AtomStamp* currentAtom;
181		BondStamp* currentBond;
182		BendStamp* currentBend;
183		TorsionStamp* currentTorsion;
184	+	RigidBodyStamp* currentRigidBody;
185
186		bond_pair* theBonds;
187		bend_set* theBends;
188		torsion_set* theTorsions;
189
190	+	set<int> skipList;
191	+
192	+	double phi, theta, psi;
193	+
194		//init the forceField paramters
195
196		the_ff->readParams();
197
187	–
198		// init the atoms
199
200	<	double phi, theta, psi;
191	<	double sux, suy, suz;
192	<	double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
193	<	double ux, uy, uz, u, uSqr;
200	>	int nMembers, nNew, rb1, rb2;
201
202		for (k = 0; k < nInfo; k++){
203		the_ff->setSimInfo(&(info[k]));
204
205		atomOffset = 0;
206	<	excludeOffset = 0;
206	>
207		for (i = 0; i < info[k].n_mol; i++){
208		stampID = info[k].molecules[i].getStampID();
209
#	Line 204 \| Line 211 \| void SimSetup::makeMolecules(void){
211		molInfo.nBonds = comp_stamps[stampID]->getNBonds();
212		molInfo.nBends = comp_stamps[stampID]->getNBends();
213		molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
214	<	molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
215	<
214	>	molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
215	>
216		molInfo.myAtoms = &(info[k].atoms[atomOffset]);
210	–	molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
211	–	molInfo.myBonds = new Bond * [molInfo.nBonds];
212	–	molInfo.myBends = new Bend * [molInfo.nBends];
213	–	molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
217
218	+	if (molInfo.nBonds > 0)
219	+	molInfo.myBonds = new (Bond *) [molInfo.nBonds];
220	+	else
221	+	molInfo.myBonds = NULL;
222	+
223	+	if (molInfo.nBends > 0)
224	+	molInfo.myBends = new (Bend *) [molInfo.nBends];
225	+	else
226	+	molInfo.myBends = NULL;
227	+
228	+	if (molInfo.nTorsions > 0)
229	+	molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
230	+	else
231	+	molInfo.myTorsions = NULL;
232	+
233		theBonds = new bond_pair[molInfo.nBonds];
234		theBends = new bend_set[molInfo.nBends];
235		theTorsions = new torsion_set[molInfo.nTorsions];
236	<
236	>
237		// make the Atoms
238
239		for (j = 0; j < molInfo.nAtoms; j++){
240		currentAtom = comp_stamps[stampID]->getAtom(j);
241	+
242		if (currentAtom->haveOrientation()){
243		dAtom = new DirectionalAtom((j + atomOffset),
244		info[k].getConfiguration());
#	Line 233 \| Line 252 \| void SimSetup::makeMolecules(void){
252		phi = currentAtom->getEulerPhi() * M_PI / 180.0;
253		theta = currentAtom->getEulerTheta() * M_PI / 180.0;
254		psi = currentAtom->getEulerPsi()* M_PI / 180.0;
236	–
237	–	Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
238	–	Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
239	–	Axz = sin(theta) * sin(psi);
240	–
241	–	Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
242	–	Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
243	–	Ayz = sin(theta) * cos(psi);
244	–
245	–	Azx = sin(phi) * sin(theta);
246	–	Azy = -cos(phi) * sin(theta);
247	–	Azz = cos(theta);
255
256	<	sux = 0.0;
257	<	suy = 0.0;
251	<	suz = 1.0;
252	<
253	<	ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
254	<	uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
255	<	uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
256	<
257	<	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
258	<
259	<	u = sqrt(uSqr);
260	<	ux = ux / u;
261	<	uy = uy / u;
262	<	uz = uz / u;
263	<
264	<	dAtom->setSUx(ux);
265	<	dAtom->setSUy(uy);
266	<	dAtom->setSUz(uz);
256	>	dAtom->setUnitFrameFromEuler(phi, theta, psi);
257	>
258		}
259		else{
260	<	molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
261	<	info[k].getConfiguration());
260	>
261	>	molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
262		}
263	+
264		molInfo.myAtoms[j]->setType(currentAtom->getType());
265
266		#ifdef IS_MPI
#	Line 284 \| Line 276 \| void SimSetup::makeMolecules(void){
276		theBonds[j].a = currentBond->getA() + atomOffset;
277		theBonds[j].b = currentBond->getB() + atomOffset;
278
279	<	exI = theBonds[j].a;
280	<	exJ = theBonds[j].b;
279	>	tempI = theBonds[j].a;
280	>	tempJ = theBonds[j].b;
281
290	–	// exclude_I must always be the smaller of the pair
291	–	if (exI > exJ){
292	–	tempEx = exI;
293	–	exI = exJ;
294	–	exJ = tempEx;
295	–	}
282		#ifdef IS_MPI
283	<	tempEx = exI;
284	<	exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
285	<	tempEx = exJ;
286	<	exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
283	>	exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
284	>	exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
285	>	#else
286	>	exI = tempI + 1;
287	>	exJ = tempJ + 1;
288	>	#endif
289
290	<	info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
303	<	#else // isn't MPI
304	<
305	<	info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
306	<	#endif //is_mpi
290	>	info[k].excludes->addPair(exI, exJ);
291		}
308	–	excludeOffset += molInfo.nBonds;
292
293		//make the bends
294		for (j = 0; j < molInfo.nBends; j++){
#	Line 355 \| Line 338 \| void SimSetup::makeMolecules(void){
338		}
339		}
340
341	<	if (!theBends[j].isGhost){
342	<	exI = theBends[j].a;
343	<	exJ = theBends[j].c;
344	<	}
345	<	else{
363	<	exI = theBends[j].a;
364	<	exJ = theBends[j].b;
365	<	}
366	<
367	<	// exclude_I must always be the smaller of the pair
368	<	if (exI > exJ){
369	<	tempEx = exI;
370	<	exI = exJ;
371	<	exJ = tempEx;
372	<	}
341	>	if (theBends[j].isGhost) {
342	>
343	>	tempI = theBends[j].a;
344	>	tempJ = theBends[j].b;
345	>
346		#ifdef IS_MPI
347	<	tempEx = exI;
348	<	exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
349	<	tempEx = exJ;
350	<	exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
347	>	exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
348	>	exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
349	>	#else
350	>	exI = tempI + 1;
351	>	exJ = tempJ + 1;
352	>	#endif
353	>	info[k].excludes->addPair(exI, exJ);
354
355	<	info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
356	<	#else // isn't MPI
357	<	info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
358	<	#endif //is_mpi
355	>	} else {
356	>
357	>	tempI = theBends[j].a;
358	>	tempJ = theBends[j].b;
359	>	tempK = theBends[j].c;
360	>
361	>	#ifdef IS_MPI
362	>	exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
363	>	exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
364	>	exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
365	>	#else
366	>	exI = tempI + 1;
367	>	exJ = tempJ + 1;
368	>	exK = tempK + 1;
369	>	#endif
370	>
371	>	info[k].excludes->addPair(exI, exK);
372	>	info[k].excludes->addPair(exI, exJ);
373	>	info[k].excludes->addPair(exJ, exK);
374	>	}
375		}
384	–	excludeOffset += molInfo.nBends;
376
377		for (j = 0; j < molInfo.nTorsions; j++){
378		currentTorsion = comp_stamps[stampID]->getTorsion(j);
#	Line 390 \| Line 381 \| void SimSetup::makeMolecules(void){
381		theTorsions[j].c = currentTorsion->getC() + atomOffset;
382		theTorsions[j].d = currentTorsion->getD() + atomOffset;
383
384	<	exI = theTorsions[j].a;
385	<	exJ = theTorsions[j].d;
384	>	tempI = theTorsions[j].a;
385	>	tempJ = theTorsions[j].b;
386	>	tempK = theTorsions[j].c;
387	>	tempL = theTorsions[j].d;
388
396	–	// exclude_I must always be the smaller of the pair
397	–	if (exI > exJ){
398	–	tempEx = exI;
399	–	exI = exJ;
400	–	exJ = tempEx;
401	–	}
389		#ifdef IS_MPI
390	<	tempEx = exI;
391	<	exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
392	<	tempEx = exJ;
393	<	exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
390	>	exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
391	>	exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
392	>	exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
393	>	exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
394	>	#else
395	>	exI = tempI + 1;
396	>	exJ = tempJ + 1;
397	>	exK = tempK + 1;
398	>	exL = tempL + 1;
399	>	#endif
400
401	<	info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
402	<	#else // isn't MPI
403	<	info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
404	<	#endif //is_mpi
401	>	info[k].excludes->addPair(exI, exJ);
402	>	info[k].excludes->addPair(exI, exK);
403	>	info[k].excludes->addPair(exI, exL);
404	>	info[k].excludes->addPair(exJ, exK);
405	>	info[k].excludes->addPair(exJ, exL);
406	>	info[k].excludes->addPair(exK, exL);
407		}
413	–	excludeOffset += molInfo.nTorsions;
408
409	+	for (j = 0; j < molInfo.nRigidBodies; j++){
410
411	<	// send the arrays off to the forceField for init.
411	>	currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
412	>	nMembers = currentRigidBody->getNMembers();
413
414	+	// Create the Rigid Body:
415	+
416	+	myRB = new RigidBody();
417	+
418	+	for (rb1 = 0; rb1 < nMembers; rb1++) {
419	+
420	+	// molI is atom numbering inside this molecule
421	+	molI = currentRigidBody->getMember(rb1);
422	+
423	+	// tempI is atom numbering on local processor
424	+	tempI = molI + atomOffset;
425	+
426	+	// currentAtom is the AtomStamp (which we need for
427	+	// rigid body reference positions)
428	+	currentAtom = comp_stamps[stampID]->getAtom(molI);
429	+
430	+	// When we add to the rigid body, add the atom itself and
431	+	// the stamp info:
432	+
433	+	myRB->addAtom(info[k].atoms[tempI], currentAtom);
434	+
435	+	// Add this atom to the Skip List for the integrators
436	+	#ifdef IS_MPI
437	+	slI = info[k].atoms[tempI]->getGlobalIndex();
438	+	#else
439	+	slI = tempI;
440	+	#endif
441	+	skipList.insert(slI);
442	+
443	+	}
444	+
445	+	for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
446	+	for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
447	+
448	+	tempI = currentRigidBody->getMember(rb1);
449	+	tempJ = currentRigidBody->getMember(rb2);
450	+
451	+	// Some explanation is required here.
452	+	// Fortran indexing starts at 1, while c indexing starts at 0
453	+	// Also, in parallel computations, the GlobalIndex is
454	+	// used for the exclude list:
455	+
456	+	#ifdef IS_MPI
457	+	exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
458	+	exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
459	+	#else
460	+	exI = tempI + 1;
461	+	exJ = tempJ + 1;
462	+	#endif
463	+
464	+	info[k].excludes->addPair(exI, exJ);
465	+
466	+	}
467	+	}
468	+	}
469	+
470	+	// send the arrays off to the forceField for init.
471	+
472		the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
473		the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
474		the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
475		the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
476		theTorsions);
477
424	–
478		info[k].molecules[i].initialize(molInfo);
479
480
#	Line 430 \| Line 483 \| void SimSetup::makeMolecules(void){
483		delete[] theBends;
484		delete[] theTorsions;
485		}
486	+
487	+	// build up the integrableObjects vector:
488	+
489	+	for (i = 0; i < info[k].n_atoms; i++) {
490	+
491	+	#ifdef IS_MPI
492	+	slI = info[k].atoms[i]->getGlobalIndex();
493	+	#else
494	+	slI = i;
495	+	#endif
496	+
497	+	if (skipList.find(slI) == skipList.end()) {
498	+	mySD = (StuntDouble *) info[k].atoms[i];
499	+	info[k].integrableObjects.push_back(mySD);
500	+	}
501	+	}
502	+	for (i = 0; i < info[k].rigidBodies.size(); i++) {
503	+	mySD = (StuntDouble *) info[k].rigidBodies[i];
504	+	info[k].integrableObjects.push_back(mySD);
505	+	}
506	+
507		}
508
509		#ifdef IS_MPI
#	Line 439 \| Line 513 \| void SimSetup::makeMolecules(void){
513
514		// clean up the forcefield
515
516	<	the_ff->calcRcut();
516	>	if (!globals->haveLJrcut()){
517	>
518	>	the_ff->calcRcut();
519	>
520	>	} else {
521	>
522	>	the_ff->setRcut( globals->getLJrcut() );
523	>	}
524	>
525		the_ff->cleanMe();
526		}
527
#	Line 1221 \| Line 1303 \| void SimSetup::calcSysValues(void){
1303		tot_bonds = 0;
1304		tot_bends = 0;
1305		tot_torsions = 0;
1306	+	tot_rigid = 0;
1307		for (i = 0; i < n_components; i++){
1308		tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1309		tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1310		tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1311		tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1312	+	tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1313		}
1314	<
1314	>
1315		tot_SRI = tot_bonds + tot_bends + tot_torsions;
1316		molMembershipArray = new int[tot_atoms];
1317
#	Line 1249 \| Line 1333 \| void SimSetup::mpiMolDivide(void){
1333		int i, j, k;
1334		int localMol, allMol;
1335		int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1336	+	int local_rigid;
1337
1338		mpiSim = new mpiSimulation(info);
1339
#	Line 1265 \| Line 1350 \| void SimSetup::mpiMolDivide(void){
1350		local_bonds = 0;
1351		local_bends = 0;
1352		local_torsions = 0;
1353	+	local_rigid = 0;
1354		globalAtomIndex = 0;
1355
1270	–
1356		for (i = 0; i < n_components; i++){
1357		for (j = 0; j < components_nmol[i]; j++){
1358		if (mol2proc[allMol] == worldRank){
#	Line 1275 \| Line 1360 \| void SimSetup::mpiMolDivide(void){
1360		local_bonds += comp_stamps[i]->getNBonds();
1361		local_bends += comp_stamps[i]->getNBends();
1362		local_torsions += comp_stamps[i]->getNTorsions();
1363	+	local_rigid += comp_stamps[i]->getNRigidBodies();
1364		localMol++;
1365		}
1366		for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
#	Line 1288 \| Line 1374 \| void SimSetup::mpiMolDivide(void){
1374		local_SRI = local_bonds + local_bends + local_torsions;
1375
1376		info[0].n_atoms = mpiSim->getMyNlocal();
1377	+
1378
1379		if (local_atoms != info[0].n_atoms){
1380		sprintf(painCave.errMsg,
#	Line 1320 \| Line 1407 \| void SimSetup::makeSysArrays(void){
1407
1408		Atom** the_atoms;
1409		Molecule* the_molecules;
1323	–	Exclude** the_excludes;
1410
1325	–
1411		for (l = 0; l < nInfo; l++){
1412		// create the atom and short range interaction arrays
1413
#	Line 1365 \| Line 1450 \| void SimSetup::makeSysArrays(void){
1450
1451		#endif // is_mpi
1452
1453	<
1454	<	if (info[l].n_SRI){
1455	<	Exclude::createArray(info[l].n_SRI);
1371	<	the_excludes = new Exclude * [info[l].n_SRI];
1372	<	for (int ex = 0; ex < info[l].n_SRI; ex++){
1373	<	the_excludes[ex] = new Exclude(ex);
1374	<	}
1375	<	info[l].globalExcludes = new int;
1376	<	info[l].n_exclude = info[l].n_SRI;
1377	<	}
1378	<	else{
1379	<	Exclude::createArray(1);
1380	<	the_excludes = new Exclude * ;
1381	<	the_excludes[0] = new Exclude(0);
1382	<	the_excludes[0]->setPair(0, 0);
1383	<	info[l].globalExcludes = new int;
1384	<	info[l].globalExcludes[0] = 0;
1385	<	info[l].n_exclude = 0;
1386	<	}
1387	<
1453	>	info[l].globalExcludes = new int;
1454	>	info[l].globalExcludes[0] = 0;
1455	>
1456		// set the arrays into the SimInfo object
1457
1458		info[l].atoms = the_atoms;
1459		info[l].molecules = the_molecules;
1460		info[l].nGlobalExcludes = 0;
1393	–	info[l].excludes = the_excludes;
1461
1462		the_ff->setSimInfo(info);
1463		}

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents): Revision 1096 by tim, Wed Mar 17 14:22:59 2004 UTC vs. Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 1096 by tim, Wed Mar 17 14:22:59 2004 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC