[ViewVC] Diff of: OpenMD/trunk/src/parallel/ForceMatrixDecomposition.cpp

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1567 by gezelter, Tue May 24 21:24:45 2011 UTC vs.
Revision 1581 by gezelter, Mon Jun 13 22:13:12 2011 UTC

#	Line 42 \| Line 42
42		#include "math/SquareMatrix3.hpp"
43		#include "nonbonded/NonBondedInteraction.hpp"
44		#include "brains/SnapshotManager.hpp"
45	+	#include "brains/PairList.hpp"
46
47		using namespace std;
48		namespace OpenMD {
#	Line 54 \| Line 55 \| namespace OpenMD {
55		void ForceMatrixDecomposition::distributeInitialData() {
56		snap_ = sman_->getCurrentSnapshot();
57		storageLayout_ = sman_->getStorageLayout();
58	+	ff_ = info_->getForceField();
59		nLocal_ = snap_->getNumberOfAtoms();
58	–	nGroups_ = snap_->getNumberOfCutoffGroups();
60
61	+	nGroups_ = info_->getNLocalCutoffGroups();
62	+	cerr << "in dId, nGroups = " << nGroups_ << "\n";
63	+	// gather the information for atomtype IDs (atids):
64	+	identsLocal = info_->getIdentArray();
65	+	AtomLocalToGlobal = info_->getGlobalAtomIndices();
66	+	cgLocalToGlobal = info_->getGlobalGroupIndices();
67	+	vector<int> globalGroupMembership = info_->getGlobalGroupMembership();
68	+	massFactors = info_->getMassFactors();
69	+	PairList excludes = info_->getExcludedInteractions();
70	+	PairList oneTwo = info_->getOneTwoInteractions();
71	+	PairList oneThree = info_->getOneThreeInteractions();
72	+	PairList oneFour = info_->getOneFourInteractions();
73	+
74		#ifdef IS_MPI
75
76		AtomCommIntRow = new Communicator<Row,int>(nLocal_);
77		AtomCommRealRow = new Communicator<Row,RealType>(nLocal_);
78		AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_);
79		AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_);
80	+	AtomCommPotRow = new Communicator<Row,potVec>(nLocal_);
81
82		AtomCommIntColumn = new Communicator<Column,int>(nLocal_);
83		AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_);
84		AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_);
85		AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_);
86	+	AtomCommPotColumn = new Communicator<Column,potVec>(nLocal_);
87
88		cgCommIntRow = new Communicator<Row,int>(nGroups_);
89		cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_);
#	Line 88 \| Line 104 \| namespace OpenMD {
104		cgRowData.setStorageLayout(DataStorage::dslPosition);
105		cgColData.resize(nGroupsInCol_);
106		cgColData.setStorageLayout(DataStorage::dslPosition);
107	+
108	+	identsRow.resize(nAtomsInRow_);
109	+	identsCol.resize(nAtomsInCol_);
110
92	–	vector<vector<RealType> > pot_row(N_INTERACTION_FAMILIES,
93	–	vector<RealType> (nAtomsInRow_, 0.0));
94	–	vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES,
95	–	vector<RealType> (nAtomsInCol_, 0.0));
96	–
97	–
98	–	vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0);
99	–
100	–	// gather the information for atomtype IDs (atids):
101	–	vector<int> identsLocal = info_->getIdentArray();
102	–	identsRow.reserve(nAtomsInRow_);
103	–	identsCol.reserve(nAtomsInCol_);
104	–
111		AtomCommIntRow->gather(identsLocal, identsRow);
112		AtomCommIntColumn->gather(identsLocal, identsCol);
113
108	–	AtomLocalToGlobal = info_->getGlobalAtomIndices();
114		AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
115		AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal);
116
112	–	cgLocalToGlobal = info_->getGlobalGroupIndices();
117		cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal);
118		cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal);
119
120	<	// still need:
121	<	// topoDist
122	<	// exclude
120	>	AtomCommRealRow->gather(massFactors, massFactorsRow);
121	>	AtomCommRealColumn->gather(massFactors, massFactorsCol);
122	>
123	>	groupListRow_.clear();
124	>	groupListRow_.resize(nGroupsInRow_);
125	>	for (int i = 0; i < nGroupsInRow_; i++) {
126	>	int gid = cgRowToGlobal[i];
127	>	for (int j = 0; j < nAtomsInRow_; j++) {
128	>	int aid = AtomRowToGlobal[j];
129	>	if (globalGroupMembership[aid] == gid)
130	>	groupListRow_[i].push_back(j);
131	>	}
132	>	}
133	>
134	>	groupListCol_.clear();
135	>	groupListCol_.resize(nGroupsInCol_);
136	>	for (int i = 0; i < nGroupsInCol_; i++) {
137	>	int gid = cgColToGlobal[i];
138	>	for (int j = 0; j < nAtomsInCol_; j++) {
139	>	int aid = AtomColToGlobal[j];
140	>	if (globalGroupMembership[aid] == gid)
141	>	groupListCol_[i].push_back(j);
142	>	}
143	>	}
144	>
145	>	skipsForAtom.clear();
146	>	skipsForAtom.resize(nAtomsInRow_);
147	>	toposForAtom.clear();
148	>	toposForAtom.resize(nAtomsInRow_);
149	>	topoDist.clear();
150	>	topoDist.resize(nAtomsInRow_);
151	>	for (int i = 0; i < nAtomsInRow_; i++) {
152	>	int iglob = AtomRowToGlobal[i];
153	>
154	>	for (int j = 0; j < nAtomsInCol_; j++) {
155	>	int jglob = AtomColToGlobal[j];
156	>
157	>	if (excludes.hasPair(iglob, jglob))
158	>	skipsForAtom[i].push_back(j);
159	>
160	>	if (oneTwo.hasPair(iglob, jglob)) {
161	>	toposForAtom[i].push_back(j);
162	>	topoDist[i].push_back(1);
163	>	} else {
164	>	if (oneThree.hasPair(iglob, jglob)) {
165	>	toposForAtom[i].push_back(j);
166	>	topoDist[i].push_back(2);
167	>	} else {
168	>	if (oneFour.hasPair(iglob, jglob)) {
169	>	toposForAtom[i].push_back(j);
170	>	topoDist[i].push_back(3);
171	>	}
172	>	}
173	>	}
174	>	}
175	>	}
176	>
177	>	#endif
178	>
179	>	groupList_.clear();
180	>	groupList_.resize(nGroups_);
181	>	for (int i = 0; i < nGroups_; i++) {
182	>	int gid = cgLocalToGlobal[i];
183	>	for (int j = 0; j < nLocal_; j++) {
184	>	int aid = AtomLocalToGlobal[j];
185	>	if (globalGroupMembership[aid] == gid) {
186	>	groupList_[i].push_back(j);
187	>	}
188	>	}
189	>	}
190	>
191	>	skipsForAtom.clear();
192	>	skipsForAtom.resize(nLocal_);
193	>	toposForAtom.clear();
194	>	toposForAtom.resize(nLocal_);
195	>	topoDist.clear();
196	>	topoDist.resize(nLocal_);
197	>
198	>	for (int i = 0; i < nLocal_; i++) {
199	>	int iglob = AtomLocalToGlobal[i];
200	>
201	>	for (int j = 0; j < nLocal_; j++) {
202	>	int jglob = AtomLocalToGlobal[j];
203	>
204	>	if (excludes.hasPair(iglob, jglob))
205	>	skipsForAtom[i].push_back(j);
206	>
207	>	if (oneTwo.hasPair(iglob, jglob)) {
208	>	toposForAtom[i].push_back(j);
209	>	topoDist[i].push_back(1);
210	>	} else {
211	>	if (oneThree.hasPair(iglob, jglob)) {
212	>	toposForAtom[i].push_back(j);
213	>	topoDist[i].push_back(2);
214	>	} else {
215	>	if (oneFour.hasPair(iglob, jglob)) {
216	>	toposForAtom[i].push_back(j);
217	>	topoDist[i].push_back(3);
218	>	}
219	>	}
220	>	}
221	>	}
222	>	}
223	>
224	>	createGtypeCutoffMap();
225	>	}
226	>
227	>	void ForceMatrixDecomposition::createGtypeCutoffMap() {
228	>
229	>	RealType tol = 1e-6;
230	>	RealType rc;
231	>	int atid;
232	>	set<AtomType*> atypes = info_->getSimulatedAtomTypes();
233	>	vector<RealType> atypeCutoff;
234	>	atypeCutoff.resize( atypes.size() );
235	>
236	>	for (set<AtomType*>::iterator at = atypes.begin();
237	>	at != atypes.end(); ++at){
238	>	rc = interactionMan_->getSuggestedCutoffRadius(*at);
239	>	atid = (*at)->getIdent();
240	>	atypeCutoff[atid] = rc;
241	>	}
242	>
243	>	vector<RealType> gTypeCutoffs;
244	>
245	>	// first we do a single loop over the cutoff groups to find the
246	>	// largest cutoff for any atypes present in this group.
247	>	#ifdef IS_MPI
248	>	vector<RealType> groupCutoffRow(nGroupsInRow_, 0.0);
249	>	groupRowToGtype.resize(nGroupsInRow_);
250	>	for (int cg1 = 0; cg1 < nGroupsInRow_; cg1++) {
251	>	vector<int> atomListRow = getAtomsInGroupRow(cg1);
252	>	for (vector<int>::iterator ia = atomListRow.begin();
253	>	ia != atomListRow.end(); ++ia) {
254	>	int atom1 = (*ia);
255	>	atid = identsRow[atom1];
256	>	if (atypeCutoff[atid] > groupCutoffRow[cg1]) {
257	>	groupCutoffRow[cg1] = atypeCutoff[atid];
258	>	}
259	>	}
260	>
261	>	bool gTypeFound = false;
262	>	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
263	>	if (abs(groupCutoffRow[cg1] - gTypeCutoffs[gt]) < tol) {
264	>	groupRowToGtype[cg1] = gt;
265	>	gTypeFound = true;
266	>	}
267	>	}
268	>	if (!gTypeFound) {
269	>	gTypeCutoffs.push_back( groupCutoffRow[cg1] );
270	>	groupRowToGtype[cg1] = gTypeCutoffs.size() - 1;
271	>	}
272	>
273	>	}
274	>	vector<RealType> groupCutoffCol(nGroupsInCol_, 0.0);
275	>	groupColToGtype.resize(nGroupsInCol_);
276	>	for (int cg2 = 0; cg2 < nGroupsInCol_; cg2++) {
277	>	vector<int> atomListCol = getAtomsInGroupColumn(cg2);
278	>	for (vector<int>::iterator jb = atomListCol.begin();
279	>	jb != atomListCol.end(); ++jb) {
280	>	int atom2 = (*jb);
281	>	atid = identsCol[atom2];
282	>	if (atypeCutoff[atid] > groupCutoffCol[cg2]) {
283	>	groupCutoffCol[cg2] = atypeCutoff[atid];
284	>	}
285	>	}
286	>	bool gTypeFound = false;
287	>	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
288	>	if (abs(groupCutoffCol[cg2] - gTypeCutoffs[gt]) < tol) {
289	>	groupColToGtype[cg2] = gt;
290	>	gTypeFound = true;
291	>	}
292	>	}
293	>	if (!gTypeFound) {
294	>	gTypeCutoffs.push_back( groupCutoffCol[cg2] );
295	>	groupColToGtype[cg2] = gTypeCutoffs.size() - 1;
296	>	}
297	>	}
298	>	#else
299	>
300	>	vector<RealType> groupCutoff(nGroups_, 0.0);
301	>	groupToGtype.resize(nGroups_);
302	>
303	>	cerr << "nGroups = " << nGroups_ << "\n";
304	>	for (int cg1 = 0; cg1 < nGroups_; cg1++) {
305	>
306	>	groupCutoff[cg1] = 0.0;
307	>	vector<int> atomList = getAtomsInGroupRow(cg1);
308	>
309	>	for (vector<int>::iterator ia = atomList.begin();
310	>	ia != atomList.end(); ++ia) {
311	>	int atom1 = (*ia);
312	>	atid = identsLocal[atom1];
313	>	if (atypeCutoff[atid] > groupCutoff[cg1]) {
314	>	groupCutoff[cg1] = atypeCutoff[atid];
315	>	}
316	>	}
317	>
318	>	bool gTypeFound = false;
319	>	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
320	>	if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) {
321	>	groupToGtype[cg1] = gt;
322	>	gTypeFound = true;
323	>	}
324	>	}
325	>	if (!gTypeFound) {
326	>	gTypeCutoffs.push_back( groupCutoff[cg1] );
327	>	groupToGtype[cg1] = gTypeCutoffs.size() - 1;
328	>	}
329	>	}
330		#endif
331	+
332	+	cerr << "gTypeCutoffs.size() = " << gTypeCutoffs.size() << "\n";
333	+	// Now we find the maximum group cutoff value present in the simulation
334	+
335	+	RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end());
336	+
337	+	#ifdef IS_MPI
338	+	MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX);
339	+	#endif
340	+
341	+	RealType tradRcut = groupMax;
342	+
343	+	for (int i = 0; i < gTypeCutoffs.size(); i++) {
344	+	for (int j = 0; j < gTypeCutoffs.size(); j++) {
345	+	RealType thisRcut;
346	+	switch(cutoffPolicy_) {
347	+	case TRADITIONAL:
348	+	thisRcut = tradRcut;
349	+	break;
350	+	case MIX:
351	+	thisRcut = 0.5 * (gTypeCutoffs[i] + gTypeCutoffs[j]);
352	+	break;
353	+	case MAX:
354	+	thisRcut = max(gTypeCutoffs[i], gTypeCutoffs[j]);
355	+	break;
356	+	default:
357	+	sprintf(painCave.errMsg,
358	+	"ForceMatrixDecomposition::createGtypeCutoffMap "
359	+	"hit an unknown cutoff policy!\n");
360	+	painCave.severity = OPENMD_ERROR;
361	+	painCave.isFatal = 1;
362	+	simError();
363	+	break;
364	+	}
365	+
366	+	pair<int,int> key = make_pair(i,j);
367	+	gTypeCutoffMap[key].first = thisRcut;
368	+
369	+	if (thisRcut > largestRcut_) largestRcut_ = thisRcut;
370	+
371	+	gTypeCutoffMap[key].second = thisRcut*thisRcut;
372	+
373	+	gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
374	+
375	+	// sanity check
376	+
377	+	if (userChoseCutoff_) {
378	+	if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) {
379	+	sprintf(painCave.errMsg,
380	+	"ForceMatrixDecomposition::createGtypeCutoffMap "
381	+	"user-specified rCut does not match computed group Cutoff\n");
382	+	painCave.severity = OPENMD_ERROR;
383	+	painCave.isFatal = 1;
384	+	simError();
385	+	}
386	+	}
387	+	}
388	+	}
389		}
390	+
391	+
392	+	groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) {
393	+	int i, j;
394	+	#ifdef IS_MPI
395	+	i = groupRowToGtype[cg1];
396	+	j = groupColToGtype[cg2];
397	+	#else
398	+	i = groupToGtype[cg1];
399	+	j = groupToGtype[cg2];
400	+	#endif
401	+	return gTypeCutoffMap[make_pair(i,j)];
402	+	}
403	+
404	+	int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) {
405	+	for (int j = 0; j < toposForAtom[atom1].size(); j++) {
406	+	if (toposForAtom[atom1][j] == atom2)
407	+	return topoDist[atom1][j];
408	+	}
409	+	return 0;
410	+	}
411	+
412	+	void ForceMatrixDecomposition::zeroWorkArrays() {
413	+
414	+	for (int j = 0; j < N_INTERACTION_FAMILIES; j++) {
415	+	longRangePot_[j] = 0.0;
416	+	}
417	+
418	+	#ifdef IS_MPI
419	+	if (storageLayout_ & DataStorage::dslForce) {
420	+	fill(atomRowData.force.begin(), atomRowData.force.end(), V3Zero);
421	+	fill(atomColData.force.begin(), atomColData.force.end(), V3Zero);
422	+	}
423	+
424	+	if (storageLayout_ & DataStorage::dslTorque) {
425	+	fill(atomRowData.torque.begin(), atomRowData.torque.end(), V3Zero);
426	+	fill(atomColData.torque.begin(), atomColData.torque.end(), V3Zero);
427	+	}
428
429	+	fill(pot_row.begin(), pot_row.end(),
430	+	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
431
432	+	fill(pot_col.begin(), pot_col.end(),
433	+	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
434	+
435	+	pot_local = Vector<RealType, N_INTERACTION_FAMILIES>(0.0);
436
437	+	if (storageLayout_ & DataStorage::dslParticlePot) {
438	+	fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0);
439	+	fill(atomColData.particlePot.begin(), atomColData.particlePot.end(), 0.0);
440	+	}
441	+
442	+	if (storageLayout_ & DataStorage::dslDensity) {
443	+	fill(atomRowData.density.begin(), atomRowData.density.end(), 0.0);
444	+	fill(atomColData.density.begin(), atomColData.density.end(), 0.0);
445	+	}
446	+
447	+	if (storageLayout_ & DataStorage::dslFunctional) {
448	+	fill(atomRowData.functional.begin(), atomRowData.functional.end(), 0.0);
449	+	fill(atomColData.functional.begin(), atomColData.functional.end(), 0.0);
450	+	}
451	+
452	+	if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
453	+	fill(atomRowData.functionalDerivative.begin(),
454	+	atomRowData.functionalDerivative.end(), 0.0);
455	+	fill(atomColData.functionalDerivative.begin(),
456	+	atomColData.functionalDerivative.end(), 0.0);
457	+	}
458	+
459	+	#else
460	+
461	+	if (storageLayout_ & DataStorage::dslParticlePot) {
462	+	fill(snap_->atomData.particlePot.begin(),
463	+	snap_->atomData.particlePot.end(), 0.0);
464	+	}
465	+
466	+	if (storageLayout_ & DataStorage::dslDensity) {
467	+	fill(snap_->atomData.density.begin(),
468	+	snap_->atomData.density.end(), 0.0);
469	+	}
470	+	if (storageLayout_ & DataStorage::dslFunctional) {
471	+	fill(snap_->atomData.functional.begin(),
472	+	snap_->atomData.functional.end(), 0.0);
473	+	}
474	+	if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
475	+	fill(snap_->atomData.functionalDerivative.begin(),
476	+	snap_->atomData.functionalDerivative.end(), 0.0);
477	+	}
478	+	#endif
479	+
480	+	}
481	+
482	+
483		void ForceMatrixDecomposition::distributeData() {
484		snap_ = sman_->getCurrentSnapshot();
485		storageLayout_ = sman_->getStorageLayout();
#	Line 156 \| Line 515 \| namespace OpenMD {
515		#endif
516		}
517
518	+	/* collects information obtained during the pre-pair loop onto local
519	+	* data structures.
520	+	*/
521		void ForceMatrixDecomposition::collectIntermediateData() {
522		snap_ = sman_->getCurrentSnapshot();
523		storageLayout_ = sman_->getStorageLayout();
#	Line 167 \| Line 529 \| namespace OpenMD {
529		snap_->atomData.density);
530
531		int n = snap_->atomData.density.size();
532	<	std::vector<RealType> rho_tmp(n, 0.0);
532	>	vector<RealType> rho_tmp(n, 0.0);
533		AtomCommRealColumn->scatter(atomColData.density, rho_tmp);
534		for (int i = 0; i < n; i++)
535		snap_->atomData.density[i] += rho_tmp[i];
536		}
537		#endif
538		}
539	<
539	>
540	>	/*
541	>	* redistributes information obtained during the pre-pair loop out to
542	>	* row and column-indexed data structures
543	>	*/
544		void ForceMatrixDecomposition::distributeIntermediateData() {
545		snap_ = sman_->getCurrentSnapshot();
546		storageLayout_ = sman_->getStorageLayout();
#	Line 232 \| Line 598 \| namespace OpenMD {
598
599		nLocal_ = snap_->getNumberOfAtoms();
600
601	<	vector<vector<RealType> > pot_temp(N_INTERACTION_FAMILIES,
602	<	vector<RealType> (nLocal_, 0.0));
601	>	vector<potVec> pot_temp(nLocal_,
602	>	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
603	>
604	>	// scatter/gather pot_row into the members of my column
605	>
606	>	AtomCommPotRow->scatter(pot_row, pot_temp);
607	>
608	>	for (int ii = 0; ii < pot_temp.size(); ii++ )
609	>	pot_local += pot_temp[ii];
610
611	<	for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
612	<	AtomCommRealRow->scatter(pot_row[i], pot_temp[i]);
613	<	for (int ii = 0; ii < pot_temp[i].size(); ii++ ) {
614	<	pot_local[i] += pot_temp[i][ii];
615	<	}
616	<	}
611	>	fill(pot_temp.begin(), pot_temp.end(),
612	>	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
613	>
614	>	AtomCommPotColumn->scatter(pot_col, pot_temp);
615	>
616	>	for (int ii = 0; ii < pot_temp.size(); ii++ )
617	>	pot_local += pot_temp[ii];
618	>
619		#endif
620		}
621
622	+	int ForceMatrixDecomposition::getNAtomsInRow() {
623	+	#ifdef IS_MPI
624	+	return nAtomsInRow_;
625	+	#else
626	+	return nLocal_;
627	+	#endif
628	+	}
629	+
630	+	/**
631	+	* returns the list of atoms belonging to this group.
632	+	*/
633	+	vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){
634	+	#ifdef IS_MPI
635	+	return groupListRow_[cg1];
636	+	#else
637	+	return groupList_[cg1];
638	+	#endif
639	+	}
640	+
641	+	vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){
642	+	#ifdef IS_MPI
643	+	return groupListCol_[cg2];
644	+	#else
645	+	return groupList_[cg2];
646	+	#endif
647	+	}
648
649		Vector3d ForceMatrixDecomposition::getIntergroupVector(int cg1, int cg2){
650		Vector3d d;
#	Line 285 \| Line 686 \| namespace OpenMD {
686		snap_->wrapVector(d);
687		return d;
688		}
689	+
690	+	RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) {
691	+	#ifdef IS_MPI
692	+	return massFactorsRow[atom1];
693	+	#else
694	+	return massFactors[atom1];
695	+	#endif
696	+	}
697	+
698	+	RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) {
699	+	#ifdef IS_MPI
700	+	return massFactorsCol[atom2];
701	+	#else
702	+	return massFactors[atom2];
703	+	#endif
704	+
705	+	}
706
707		Vector3d ForceMatrixDecomposition::getInteratomicVector(int atom1, int atom2){
708		Vector3d d;
#	Line 299 \| Line 717 \| namespace OpenMD {
717		return d;
718		}
719
720	+	vector<int> ForceMatrixDecomposition::getSkipsForAtom(int atom1) {
721	+	return skipsForAtom[atom1];
722	+	}
723	+
724	+	/**
725	+	* There are a number of reasons to skip a pair or a
726	+	* particle. Mostly we do this to exclude atoms who are involved in
727	+	* short range interactions (bonds, bends, torsions), but we also
728	+	* need to exclude some overcounted interactions that result from
729	+	* the parallel decomposition.
730	+	*/
731	+	bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
732	+	int unique_id_1, unique_id_2;
733	+
734	+	#ifdef IS_MPI
735	+	// in MPI, we have to look up the unique IDs for each atom
736	+	unique_id_1 = AtomRowToGlobal[atom1];
737	+	unique_id_2 = AtomColToGlobal[atom2];
738	+
739	+	// this situation should only arise in MPI simulations
740	+	if (unique_id_1 == unique_id_2) return true;
741	+
742	+	// this prevents us from doing the pair on multiple processors
743	+	if (unique_id_1 < unique_id_2) {
744	+	if ((unique_id_1 + unique_id_2) % 2 == 0) return true;
745	+	} else {
746	+	if ((unique_id_1 + unique_id_2) % 2 == 1) return true;
747	+	}
748	+	#else
749	+	// in the normal loop, the atom numbers are unique
750	+	unique_id_1 = atom1;
751	+	unique_id_2 = atom2;
752	+	#endif
753	+
754	+	for (vector<int>::iterator i = skipsForAtom[atom1].begin();
755	+	i != skipsForAtom[atom1].end(); ++i) {
756	+	if ( (*i) == unique_id_2 ) return true;
757	+	}
758	+
759	+	}
760	+
761	+
762		void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){
763		#ifdef IS_MPI
764		atomRowData.force[atom1] += fg;
#	Line 316 \| Line 776 \| namespace OpenMD {
776		}
777
778		// filling interaction blocks with pointers
779	<	InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) {
780	<	InteractionData idat;
321	<
779	>	void ForceMatrixDecomposition::fillInteractionData(InteractionData idat,
780	>	int atom1, int atom2) {
781		#ifdef IS_MPI
782	+
783	+	idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
784	+	ff_->getAtomType(identsCol[atom2]) );
785	+
786		if (storageLayout_ & DataStorage::dslAmat) {
787		idat.A1 = &(atomRowData.aMat[atom1]);
788		idat.A2 = &(atomColData.aMat[atom2]);
#	Line 340 \| Line 803 \| namespace OpenMD {
803		idat.rho2 = &(atomColData.density[atom2]);
804		}
805
806	+	if (storageLayout_ & DataStorage::dslFunctional) {
807	+	idat.frho1 = &(atomRowData.functional[atom1]);
808	+	idat.frho2 = &(atomColData.functional[atom2]);
809	+	}
810	+
811		if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
812		idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
813		idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
814		}
815	+
816	+	if (storageLayout_ & DataStorage::dslParticlePot) {
817	+	idat.particlePot1 = &(atomRowData.particlePot[atom1]);
818	+	idat.particlePot2 = &(atomColData.particlePot[atom2]);
819	+	}
820	+
821		#else
822	+
823	+	idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]),
824	+	ff_->getAtomType(identsLocal[atom2]) );
825	+
826		if (storageLayout_ & DataStorage::dslAmat) {
827		idat.A1 = &(snap_->atomData.aMat[atom1]);
828		idat.A2 = &(snap_->atomData.aMat[atom2]);
#	Line 365 \| Line 843 \| namespace OpenMD {
843		idat.rho2 = &(snap_->atomData.density[atom2]);
844		}
845
846	+	if (storageLayout_ & DataStorage::dslFunctional) {
847	+	idat.frho1 = &(snap_->atomData.functional[atom1]);
848	+	idat.frho2 = &(snap_->atomData.functional[atom2]);
849	+	}
850	+
851		if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
852		idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
853		idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
854		}
855	+
856	+	if (storageLayout_ & DataStorage::dslParticlePot) {
857	+	idat.particlePot1 = &(snap_->atomData.particlePot[atom1]);
858	+	idat.particlePot2 = &(snap_->atomData.particlePot[atom2]);
859	+	}
860	+
861		#endif
373	–	return idat;
862		}
863
864	<	InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){
864	>
865	>	void ForceMatrixDecomposition::unpackInteractionData(InteractionData idat, int atom1, int atom2) {
866	>	#ifdef IS_MPI
867	>	pot_row[atom1] += 0.5 * *(idat.pot);
868	>	pot_col[atom2] += 0.5 * *(idat.pot);
869
870	<	InteractionData idat;
870	>	atomRowData.force[atom1] += *(idat.f1);
871	>	atomColData.force[atom2] -= *(idat.f1);
872	>	#else
873	>	longRangePot_ += *(idat.pot);
874	>
875	>	snap_->atomData.force[atom1] += *(idat.f1);
876	>	snap_->atomData.force[atom2] -= *(idat.f1);
877	>	#endif
878	>
879	>	}
880	>
881	>
882	>	void ForceMatrixDecomposition::fillSkipData(InteractionData idat,
883	>	int atom1, int atom2) {
884		#ifdef IS_MPI
885	+	idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
886	+	ff_->getAtomType(identsCol[atom2]) );
887	+
888		if (storageLayout_ & DataStorage::dslElectroFrame) {
889		idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
890		idat.eFrame2 = &(atomColData.electroFrame[atom2]);
#	Line 385 \| Line 893 \| namespace OpenMD {
893		idat.t1 = &(atomRowData.torque[atom1]);
894		idat.t2 = &(atomColData.torque[atom2]);
895		}
388	–	if (storageLayout_ & DataStorage::dslForce) {
389	–	idat.t1 = &(atomRowData.force[atom1]);
390	–	idat.t2 = &(atomColData.force[atom2]);
391	–	}
896		#else
897	+	idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]),
898	+	ff_->getAtomType(identsLocal[atom2]) );
899	+
900		if (storageLayout_ & DataStorage::dslElectroFrame) {
901		idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
902		idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
#	Line 398 \| Line 905 \| namespace OpenMD {
905		idat.t1 = &(snap_->atomData.torque[atom1]);
906		idat.t2 = &(snap_->atomData.torque[atom2]);
907		}
908	<	if (storageLayout_ & DataStorage::dslForce) {
402	<	idat.t1 = &(snap_->atomData.force[atom1]);
403	<	idat.t2 = &(snap_->atomData.force[atom2]);
404	<	}
405	<	#endif
406	<
407	<	}
408	<
409	<	SelfData ForceMatrixDecomposition::fillSelfData(int atom1) {
410	<	SelfData sdat;
411	<	// Still Missing atype, skippedCharge, potVec pot,
412	<	if (storageLayout_ & DataStorage::dslElectroFrame) {
413	<	sdat.eFrame = &(snap_->atomData.electroFrame[atom1]);
414	<	}
415	<
416	<	if (storageLayout_ & DataStorage::dslTorque) {
417	<	sdat.t = &(snap_->atomData.torque[atom1]);
418	<	}
419	<
420	<	if (storageLayout_ & DataStorage::dslDensity) {
421	<	sdat.rho = &(snap_->atomData.density[atom1]);
422	<	}
423	<
424	<	if (storageLayout_ & DataStorage::dslFunctional) {
425	<	sdat.frho = &(snap_->atomData.functional[atom1]);
426	<	}
427	<
428	<	if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
429	<	sdat.dfrhodrho = &(snap_->atomData.functionalDerivative[atom1]);
430	<	}
431	<
432	<	return sdat;
908	>	#endif
909		}
910
435	–
436	–
911		/*
912		* buildNeighborList
913		*
#	Line 443 \| Line 917 \| namespace OpenMD {
917		vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
918
919		vector<pair<int, int> > neighborList;
920	+	groupCutoffs cuts;
921		#ifdef IS_MPI
922	<	CellListRow.clear();
923	<	CellListCol.clear();
922	>	cellListRow_.clear();
923	>	cellListCol_.clear();
924		#else
925	<	CellList.clear();
925	>	cellList_.clear();
926		#endif
927
928	<	// dangerous to not do error checking.
454	<	RealType skinThickness_ = info_->getSimParams()->getSkinThickness();
455	<	RealType rCut_;
456	<
457	<	RealType rList_ = (rCut_ + skinThickness_);
928	>	RealType rList_ = (largestRcut_ + skinThickness_);
929		RealType rl2 = rList_ * rList_;
930		Snapshot* snap_ = sman_->getCurrentSnapshot();
931		Mat3x3d Hmat = snap_->getHmat();
932		Vector3d Hx = Hmat.getColumn(0);
933		Vector3d Hy = Hmat.getColumn(1);
934		Vector3d Hz = Hmat.getColumn(2);
464	–	Vector3i nCells;
935
936	<	nCells.x() = (int) ( Hx.length() )/ rList_;
937	<	nCells.y() = (int) ( Hy.length() )/ rList_;
938	<	nCells.z() = (int) ( Hz.length() )/ rList_;
936	>	nCells_.x() = (int) ( Hx.length() )/ rList_;
937	>	nCells_.y() = (int) ( Hy.length() )/ rList_;
938	>	nCells_.z() = (int) ( Hz.length() )/ rList_;
939
940		Mat3x3d invHmat = snap_->getInvHmat();
941		Vector3d rs, scaled, dr;
942		Vector3i whichCell;
943		int cellIndex;
944	+	int nCtot = nCells_.x() * nCells_.y() * nCells_.z();
945
946	+	cerr << "flag1\n";
947		#ifdef IS_MPI
948	+	cellListRow_.resize(nCtot);
949	+	cellListCol_.resize(nCtot);
950	+	#else
951	+	cellList_.resize(nCtot);
952	+	#endif
953	+	cerr << "flag2\n";
954	+	#ifdef IS_MPI
955		for (int i = 0; i < nGroupsInRow_; i++) {
956		rs = cgRowData.position[i];
957	+
958		// scaled positions relative to the box vectors
959		scaled = invHmat * rs;
960	+
961		// wrap the vector back into the unit box by subtracting integer box
962		// numbers
963	<	for (int j = 0; j < 3; j++)
963	>	for (int j = 0; j < 3; j++) {
964		scaled[j] -= roundMe(scaled[j]);
965	+	scaled[j] += 0.5;
966	+	}
967
968		// find xyz-indices of cell that cutoffGroup is in.
969	<	whichCell.x() = nCells.x() * scaled.x();
970	<	whichCell.y() = nCells.y() * scaled.y();
971	<	whichCell.z() = nCells.z() * scaled.z();
969	>	whichCell.x() = nCells_.x() * scaled.x();
970	>	whichCell.y() = nCells_.y() * scaled.y();
971	>	whichCell.z() = nCells_.z() * scaled.z();
972
973		// find single index of this cell:
974	<	cellIndex = Vlinear(whichCell, nCells);
974	>	cellIndex = Vlinear(whichCell, nCells_);
975	>
976		// add this cutoff group to the list of groups in this cell;
977	<	CellListRow[cellIndex].push_back(i);
977	>	cellListRow_[cellIndex].push_back(i);
978		}
979
980		for (int i = 0; i < nGroupsInCol_; i++) {
981		rs = cgColData.position[i];
982	+
983		// scaled positions relative to the box vectors
984		scaled = invHmat * rs;
985	+
986		// wrap the vector back into the unit box by subtracting integer box
987		// numbers
988	<	for (int j = 0; j < 3; j++)
988	>	for (int j = 0; j < 3; j++) {
989		scaled[j] -= roundMe(scaled[j]);
990	+	scaled[j] += 0.5;
991	+	}
992
993		// find xyz-indices of cell that cutoffGroup is in.
994	<	whichCell.x() = nCells.x() * scaled.x();
995	<	whichCell.y() = nCells.y() * scaled.y();
996	<	whichCell.z() = nCells.z() * scaled.z();
994	>	whichCell.x() = nCells_.x() * scaled.x();
995	>	whichCell.y() = nCells_.y() * scaled.y();
996	>	whichCell.z() = nCells_.z() * scaled.z();
997
998		// find single index of this cell:
999	<	cellIndex = Vlinear(whichCell, nCells);
999	>	cellIndex = Vlinear(whichCell, nCells_);
1000	>
1001		// add this cutoff group to the list of groups in this cell;
1002	<	CellListCol[cellIndex].push_back(i);
1002	>	cellListCol_[cellIndex].push_back(i);
1003		}
1004		#else
1005		for (int i = 0; i < nGroups_; i++) {
1006		rs = snap_->cgData.position[i];
1007	+
1008		// scaled positions relative to the box vectors
1009		scaled = invHmat * rs;
1010	+
1011		// wrap the vector back into the unit box by subtracting integer box
1012		// numbers
1013	<	for (int j = 0; j < 3; j++)
1013	>	for (int j = 0; j < 3; j++) {
1014		scaled[j] -= roundMe(scaled[j]);
1015	+	scaled[j] += 0.5;
1016	+	}
1017
1018		// find xyz-indices of cell that cutoffGroup is in.
1019	<	whichCell.x() = nCells.x() * scaled.x();
1020	<	whichCell.y() = nCells.y() * scaled.y();
1021	<	whichCell.z() = nCells.z() * scaled.z();
1019	>	whichCell.x() = nCells_.x() * scaled.x();
1020	>	whichCell.y() = nCells_.y() * scaled.y();
1021	>	whichCell.z() = nCells_.z() * scaled.z();
1022
1023		// find single index of this cell:
1024	<	cellIndex = Vlinear(whichCell, nCells);
1024	>	cellIndex = Vlinear(whichCell, nCells_);
1025	>
1026		// add this cutoff group to the list of groups in this cell;
1027	<	CellList[cellIndex].push_back(i);
1027	>	cellList_[cellIndex].push_back(i);
1028		}
1029		#endif
1030
1031	<
1032	<
1033	<	for (int m1z = 0; m1z < nCells.z(); m1z++) {
540	<	for (int m1y = 0; m1y < nCells.y(); m1y++) {
541	<	for (int m1x = 0; m1x < nCells.x(); m1x++) {
1031	>	for (int m1z = 0; m1z < nCells_.z(); m1z++) {
1032	>	for (int m1y = 0; m1y < nCells_.y(); m1y++) {
1033	>	for (int m1x = 0; m1x < nCells_.x(); m1x++) {
1034		Vector3i m1v(m1x, m1y, m1z);
1035	<	int m1 = Vlinear(m1v, nCells);
544	<	for (int offset = 0; offset < nOffset_; offset++) {
545	<	Vector3i m2v = m1v + cellOffsets_[offset];
1035	>	int m1 = Vlinear(m1v, nCells_);
1036
1037	<	if (m2v.x() >= nCells.x()) {
1037	>	for (vector<Vector3i>::iterator os = cellOffsets_.begin();
1038	>	os != cellOffsets_.end(); ++os) {
1039	>
1040	>	Vector3i m2v = m1v + (*os);
1041	>
1042	>	if (m2v.x() >= nCells_.x()) {
1043		m2v.x() = 0;
1044		} else if (m2v.x() < 0) {
1045	<	m2v.x() = nCells.x() - 1;
1045	>	m2v.x() = nCells_.x() - 1;
1046		}
1047	<
1048	<	if (m2v.y() >= nCells.y()) {
1047	>
1048	>	if (m2v.y() >= nCells_.y()) {
1049		m2v.y() = 0;
1050		} else if (m2v.y() < 0) {
1051	<	m2v.y() = nCells.y() - 1;
1051	>	m2v.y() = nCells_.y() - 1;
1052		}
1053	<
1054	<	if (m2v.z() >= nCells.z()) {
1053	>
1054	>	if (m2v.z() >= nCells_.z()) {
1055		m2v.z() = 0;
1056		} else if (m2v.z() < 0) {
1057	<	m2v.z() = nCells.z() - 1;
1057	>	m2v.z() = nCells_.z() - 1;
1058		}
1059	+
1060	+	int m2 = Vlinear (m2v, nCells_);
1061
565	–	int m2 = Vlinear (m2v, nCells);
566	–
1062		#ifdef IS_MPI
1063	<	for (vector<int>::iterator j1 = CellListRow[m1].begin();
1064	<	j1 != CellListRow[m1].end(); ++j1) {
1065	<	for (vector<int>::iterator j2 = CellListCol[m2].begin();
1066	<	j2 != CellListCol[m2].end(); ++j2) {
1063	>	for (vector<int>::iterator j1 = cellListRow_[m1].begin();
1064	>	j1 != cellListRow_[m1].end(); ++j1) {
1065	>	for (vector<int>::iterator j2 = cellListCol_[m2].begin();
1066	>	j2 != cellListCol_[m2].end(); ++j2) {
1067
1068		// Always do this if we're in different cells or if
1069		// we're in the same cell and the global index of the
#	Line 577 \| Line 1072 \| namespace OpenMD {
1072		if (m2 != m1 \|\| cgColToGlobal[(j2)] < cgRowToGlobal[(j1)]) {
1073		dr = cgColData.position[(j2)] - cgRowData.position[(j1)];
1074		snap_->wrapVector(dr);
1075	<	if (dr.lengthSquare() < rl2) {
1075	>	cuts = getGroupCutoffs( (j1), (j2) );
1076	>	if (dr.lengthSquare() < cuts.third) {
1077		neighborList.push_back(make_pair((j1), (j2)));
1078		}
1079		}
1080		}
1081		}
1082		#else
1083	<	for (vector<int>::iterator j1 = CellList[m1].begin();
1084	<	j1 != CellList[m1].end(); ++j1) {
1085	<	for (vector<int>::iterator j2 = CellList[m2].begin();
1086	<	j2 != CellList[m2].end(); ++j2) {
1087	<
1083	>
1084	>	for (vector<int>::iterator j1 = cellList_[m1].begin();
1085	>	j1 != cellList_[m1].end(); ++j1) {
1086	>	for (vector<int>::iterator j2 = cellList_[m2].begin();
1087	>	j2 != cellList_[m2].end(); ++j2) {
1088	>
1089		// Always do this if we're in different cells or if
1090		// we're in the same cell and the global index of the
1091		// j2 cutoff group is less than the j1 cutoff group
#	Line 596 \| Line 1093 \| namespace OpenMD {
1093		if (m2 != m1 \|\| (j2) < (j1)) {
1094		dr = snap_->cgData.position[(j2)] - snap_->cgData.position[(j1)];
1095		snap_->wrapVector(dr);
1096	<	if (dr.lengthSquare() < rl2) {
1096	>	cuts = getGroupCutoffs( (j1), (j2) );
1097	>	if (dr.lengthSquare() < cuts.third) {
1098		neighborList.push_back(make_pair((j1), (j2)));
1099		}
1100		}
#	Line 607 \| Line 1105 \| namespace OpenMD {
1105		}
1106		}
1107		}
1108	+
1109	+	// save the local cutoff group positions for the check that is
1110	+	// done on each loop:
1111	+	saved_CG_positions_.clear();
1112	+	for (int i = 0; i < nGroups_; i++)
1113	+	saved_CG_positions_.push_back(snap_->cgData.position[i]);
1114	+
1115		return neighborList;
1116		}
1117		} //end namespace OpenMD

Diff Legend

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

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents): Revision 1567 by gezelter, Tue May 24 21:24:45 2011 UTC vs. Revision 1581 by gezelter, Mon Jun 13 22:13:12 2011 UTC

Diff Legend

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1567 by gezelter, Tue May 24 21:24:45 2011 UTC vs.
Revision 1581 by gezelter, Mon Jun 13 22:13:12 2011 UTC