ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/parallel/ForceMatrixDecomposition.cpp
(Generate patch)

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1588 by gezelter, Sat Jul 9 15:05:59 2011 UTC vs.
Revision 1592 by gezelter, Tue Jul 12 20:33:14 2011 UTC

# Line 112 | Line 112 | namespace OpenMD {
112      AtomCommIntRow->gather(idents, identsRow);
113      AtomCommIntColumn->gather(idents, identsCol);
114      
115 <    vector<int>::iterator it;
116 <    for (it = AtomLocalToGlobal.begin(); it != AtomLocalToGlobal.end(); ++it) {
117 <      cerr << "my AtomLocalToGlobal = " << (*it) << "\n";
118 <    }
115 >    // allocate memory for the parallel objects
116 >    atypesRow.resize(nAtomsInRow_);
117 >    atypesCol.resize(nAtomsInCol_);
118 >
119 >    for (int i = 0; i < nAtomsInRow_; i++)
120 >      atypesRow[i] = ff_->getAtomType(identsRow[i]);
121 >    for (int i = 0; i < nAtomsInCol_; i++)
122 >      atypesCol[i] = ff_->getAtomType(identsCol[i]);        
123 >
124 >    pot_row.resize(nAtomsInRow_);
125 >    pot_col.resize(nAtomsInCol_);
126 >
127 >    AtomRowToGlobal.resize(nAtomsInRow_);
128 >    AtomColToGlobal.resize(nAtomsInCol_);
129      AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
130      AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal);
131      
132 +    cgRowToGlobal.resize(nGroupsInRow_);
133 +    cgColToGlobal.resize(nGroupsInCol_);
134      cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal);
135      cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal);
136  
137 +    massFactorsRow.resize(nAtomsInRow_);
138 +    massFactorsCol.resize(nAtomsInCol_);
139      AtomCommRealRow->gather(massFactors, massFactorsRow);
140      AtomCommRealColumn->gather(massFactors, massFactorsCol);
141  
# Line 180 | Line 194 | namespace OpenMD {
194      }
195  
196   #endif
197 +
198 +    // allocate memory for the parallel objects
199 +    atypesLocal.resize(nLocal_);
200  
201 +    for (int i = 0; i < nLocal_; i++)
202 +      atypesLocal[i] = ff_->getAtomType(idents[i]);
203 +
204      groupList_.clear();
205      groupList_.resize(nGroups_);
206      for (int i = 0; i < nGroups_; i++) {
# Line 233 | Line 253 | namespace OpenMD {
253    void ForceMatrixDecomposition::createGtypeCutoffMap() {
254      
255      RealType tol = 1e-6;
256 +    largestRcut_ = 0.0;
257      RealType rc;
258      int atid;
259      set<AtomType*> atypes = info_->getSimulatedAtomTypes();
260 +    
261      map<int, RealType> atypeCutoff;
262        
263      for (set<AtomType*>::iterator at = atypes.begin();
# Line 243 | Line 265 | namespace OpenMD {
265        atid = (*at)->getIdent();
266        if (userChoseCutoff_)
267          atypeCutoff[atid] = userCutoff_;
268 <      else
268 >      else
269          atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at);
270      }
271 <
271 >    
272      vector<RealType> gTypeCutoffs;
273      // first we do a single loop over the cutoff groups to find the
274      // largest cutoff for any atypes present in this group.
# Line 306 | Line 328 | namespace OpenMD {
328      vector<RealType> groupCutoff(nGroups_, 0.0);
329      groupToGtype.resize(nGroups_);
330      for (int cg1 = 0; cg1 < nGroups_; cg1++) {
309 –
331        groupCutoff[cg1] = 0.0;
332        vector<int> atomList = getAtomsInGroupRow(cg1);
312 –
333        for (vector<int>::iterator ia = atomList.begin();
334             ia != atomList.end(); ++ia) {            
335          int atom1 = (*ia);
336          atid = idents[atom1];
337 <        if (atypeCutoff[atid] > groupCutoff[cg1]) {
337 >        if (atypeCutoff[atid] > groupCutoff[cg1])
338            groupCutoff[cg1] = atypeCutoff[atid];
319 –        }
339        }
340 <
340 >      
341        bool gTypeFound = false;
342        for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
343          if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) {
# Line 326 | Line 345 | namespace OpenMD {
345            gTypeFound = true;
346          }
347        }
348 <      if (!gTypeFound) {
348 >      if (!gTypeFound) {      
349          gTypeCutoffs.push_back( groupCutoff[cg1] );
350          groupToGtype[cg1] = gTypeCutoffs.size() - 1;
351        }      
# Line 335 | Line 354 | namespace OpenMD {
354  
355      // Now we find the maximum group cutoff value present in the simulation
356  
357 <    RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end());
357 >    RealType groupMax = *max_element(gTypeCutoffs.begin(),
358 >                                     gTypeCutoffs.end());
359  
360   #ifdef IS_MPI
361 <    MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX);
361 >    MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE,
362 >                              MPI::MAX);
363   #endif
364      
365      RealType tradRcut = groupMax;
# Line 368 | Line 389 | namespace OpenMD {
389  
390          pair<int,int> key = make_pair(i,j);
391          gTypeCutoffMap[key].first = thisRcut;
371 –
392          if (thisRcut > largestRcut_) largestRcut_ = thisRcut;
373 –
393          gTypeCutoffMap[key].second = thisRcut*thisRcut;
375 –        
394          gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
377 –
395          // sanity check
396          
397          if (userChoseCutoff_) {
# Line 434 | Line 451 | namespace OpenMD {
451           Vector<RealType, N_INTERACTION_FAMILIES> (0.0));  
452  
453      if (storageLayout_ & DataStorage::dslParticlePot) {    
454 <      fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0);
455 <      fill(atomColData.particlePot.begin(), atomColData.particlePot.end(), 0.0);
454 >      fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(),
455 >           0.0);
456 >      fill(atomColData.particlePot.begin(), atomColData.particlePot.end(),
457 >           0.0);
458      }
459  
460      if (storageLayout_ & DataStorage::dslDensity) {      
# Line 444 | Line 463 | namespace OpenMD {
463      }
464  
465      if (storageLayout_ & DataStorage::dslFunctional) {  
466 <      fill(atomRowData.functional.begin(), atomRowData.functional.end(), 0.0);
467 <      fill(atomColData.functional.begin(), atomColData.functional.end(), 0.0);
466 >      fill(atomRowData.functional.begin(), atomRowData.functional.end(),
467 >           0.0);
468 >      fill(atomColData.functional.begin(), atomColData.functional.end(),
469 >           0.0);
470      }
471  
472      if (storageLayout_ & DataStorage::dslFunctionalDerivative) {      
# Line 462 | Line 483 | namespace OpenMD {
483             atomColData.skippedCharge.end(), 0.0);
484      }
485  
486 < #else
487 <    
486 > #endif
487 >    // even in parallel, we need to zero out the local arrays:
488 >
489      if (storageLayout_ & DataStorage::dslParticlePot) {      
490        fill(snap_->atomData.particlePot.begin(),
491             snap_->atomData.particlePot.end(), 0.0);
# Line 485 | Line 507 | namespace OpenMD {
507        fill(snap_->atomData.skippedCharge.begin(),
508             snap_->atomData.skippedCharge.end(), 0.0);
509      }
488 – #endif
510      
511    }
512  
# Line 522 | Line 543 | namespace OpenMD {
543        AtomCommMatrixColumn->gather(snap_->atomData.electroFrame,
544                                     atomColData.electroFrame);
545      }
546 +
547   #endif      
548    }
549    
# Line 588 | Line 610 | namespace OpenMD {
610      AtomCommVectorColumn->scatter(atomColData.force, frc_tmp);
611      for (int i = 0; i < n; i++)
612        snap_->atomData.force[i] += frc_tmp[i];
613 <    
592 <    
613 >        
614      if (storageLayout_ & DataStorage::dslTorque) {
615  
616        int nt = snap_->atomData.torque.size();
# Line 613 | Line 634 | namespace OpenMD {
634  
635        AtomCommRealRow->scatter(atomRowData.skippedCharge, skch_tmp);
636        for (int i = 0; i < ns; i++) {
637 <        snap_->atomData.skippedCharge[i] = skch_tmp[i];
637 >        snap_->atomData.skippedCharge[i] += skch_tmp[i];
638          skch_tmp[i] = 0.0;
639        }
640        
# Line 824 | Line 845 | namespace OpenMD {
845      idat.excluded = excludeAtomPair(atom1, atom2);
846    
847   #ifdef IS_MPI
848 <    
849 <    idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
850 <                             ff_->getAtomType(identsCol[atom2]) );
848 >    idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]);
849 >    //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
850 >    //                         ff_->getAtomType(identsCol[atom2]) );
851      
852      if (storageLayout_ & DataStorage::dslAmat) {
853        idat.A1 = &(atomRowData.aMat[atom1]);
# Line 870 | Line 891 | namespace OpenMD {
891  
892   #else
893  
894 <    idat.atypes = make_pair( ff_->getAtomType(idents[atom1]),
895 <                             ff_->getAtomType(idents[atom2]) );
894 >    idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]);
895 >    //idat.atypes = make_pair( ff_->getAtomType(idents[atom1]),
896 >    //                         ff_->getAtomType(idents[atom2]) );
897  
898      if (storageLayout_ & DataStorage::dslAmat) {
899        idat.A1 = &(snap_->atomData.aMat[atom1]);

Diff Legend

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