ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/ForceManager.cpp
(Generate patch)

Comparing branches/development/src/brains/ForceManager.cpp (file contents):
Revision 1545 by gezelter, Fri Apr 8 21:25:19 2011 UTC vs.
Revision 1576 by gezelter, Wed Jun 8 16:05:07 2011 UTC

# Line 47 | Line 47
47   * @version 1.0
48   */
49  
50 +
51   #include "brains/ForceManager.hpp"
52   #include "primitives/Molecule.hpp"
52 – #include "UseTheForce/doForces_interface.h"
53   #define __OPENMD_C
54 – #include "UseTheForce/DarkSide/fInteractionMap.h"
54   #include "utils/simError.h"
55   #include "primitives/Bond.hpp"
56   #include "primitives/Bend.hpp"
57   #include "primitives/Torsion.hpp"
58   #include "primitives/Inversion.hpp"
59 < #include "parallel/ForceDecomposition.hpp"
60 < //#include "parallel/SerialDecomposition.hpp"
59 > #include "nonbonded/NonBondedInteraction.hpp"
60 > #include "parallel/ForceMatrixDecomposition.hpp"
61  
62   using namespace std;
63   namespace OpenMD {
64    
65    ForceManager::ForceManager(SimInfo * info) : info_(info) {
66 <
67 < #ifdef IS_MPI
69 <    decomp_ = new ForceDecomposition(info_);
70 < #else
71 <    // decomp_ = new SerialDecomposition(info);
72 < #endif
66 >    forceField_ = info_->getForceField();
67 >    fDecomp_ = new ForceMatrixDecomposition(info_);
68    }
69 +
70 +  /**
71 +   * setupCutoffs
72 +   *
73 +   * Sets the values of cutoffRadius, cutoffMethod, and cutoffPolicy
74 +   *
75 +   * cutoffRadius : realType
76 +   *  If the cutoffRadius was explicitly set, use that value.
77 +   *  If the cutoffRadius was not explicitly set:
78 +   *      Are there electrostatic atoms?  Use 12.0 Angstroms.
79 +   *      No electrostatic atoms?  Poll the atom types present in the
80 +   *      simulation for suggested cutoff values (e.g. 2.5 * sigma).
81 +   *      Use the maximum suggested value that was found.
82 +   *
83 +   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
84 +   *      If cutoffMethod was explicitly set, use that choice.
85 +   *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
86 +   *
87 +   * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
88 +   *      If cutoffPolicy was explicitly set, use that choice.
89 +   *      If cutoffPolicy was not explicitly set, use TRADITIONAL
90 +   */
91 +  void ForceManager::setupCutoffs() {
92 +    
93 +    Globals* simParams_ = info_->getSimParams();
94 +    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
95 +    
96 +    if (simParams_->haveCutoffRadius()) {
97 +      rCut_ = simParams_->getCutoffRadius();
98 +    } else {      
99 +      if (info_->usesElectrostaticAtoms()) {
100 +        sprintf(painCave.errMsg,
101 +                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
102 +                "\tOpenMD will use a default value of 12.0 angstroms"
103 +                "\tfor the cutoffRadius.\n");
104 +        painCave.isFatal = 0;
105 +        painCave.severity = OPENMD_INFO;
106 +        simError();
107 +        rCut_ = 12.0;
108 +      } else {
109 +        RealType thisCut;
110 +        set<AtomType*>::iterator i;
111 +        set<AtomType*> atomTypes;
112 +        atomTypes = info_->getSimulatedAtomTypes();        
113 +        for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
114 +          thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
115 +          rCut_ = max(thisCut, rCut_);
116 +        }
117 +        sprintf(painCave.errMsg,
118 +                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
119 +                "\tOpenMD will use %lf angstroms.\n",
120 +                rCut_);
121 +        painCave.isFatal = 0;
122 +        painCave.severity = OPENMD_INFO;
123 +        simError();
124 +      }            
125 +    }
126 +
127 +    map<string, CutoffMethod> stringToCutoffMethod;
128 +    stringToCutoffMethod["HARD"] = HARD;
129 +    stringToCutoffMethod["SWITCHED"] = SWITCHED;
130 +    stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;    
131 +    stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
132    
133 <  void ForceManager::calcForces() {
133 >    if (simParams_->haveCutoffMethod()) {
134 >      string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
135 >      map<string, CutoffMethod>::iterator i;
136 >      i = stringToCutoffMethod.find(cutMeth);
137 >      if (i == stringToCutoffMethod.end()) {
138 >        sprintf(painCave.errMsg,
139 >                "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
140 >                "\tShould be one of: "
141 >                "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
142 >                cutMeth.c_str());
143 >        painCave.isFatal = 1;
144 >        painCave.severity = OPENMD_ERROR;
145 >        simError();
146 >      } else {
147 >        cutoffMethod_ = i->second;
148 >      }
149 >    } else {
150 >      sprintf(painCave.errMsg,
151 >              "ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n"
152 >              "\tOpenMD will use SHIFTED_FORCE.\n");
153 >      painCave.isFatal = 0;
154 >      painCave.severity = OPENMD_INFO;
155 >      simError();
156 >      cutoffMethod_ = SHIFTED_FORCE;        
157 >    }
158 >
159 >    map<string, CutoffPolicy> stringToCutoffPolicy;
160 >    stringToCutoffPolicy["MIX"] = MIX;
161 >    stringToCutoffPolicy["MAX"] = MAX;
162 >    stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;    
163 >
164 >    std::string cutPolicy;
165 >    if (forceFieldOptions_.haveCutoffPolicy()){
166 >      cutPolicy = forceFieldOptions_.getCutoffPolicy();
167 >    }else if (simParams_->haveCutoffPolicy()) {
168 >      cutPolicy = simParams_->getCutoffPolicy();
169 >    }
170 >
171 >    if (!cutPolicy.empty()){
172 >      toUpper(cutPolicy);
173 >      map<string, CutoffPolicy>::iterator i;
174 >      i = stringToCutoffPolicy.find(cutPolicy);
175 >
176 >      if (i == stringToCutoffPolicy.end()) {
177 >        sprintf(painCave.errMsg,
178 >                "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
179 >                "\tShould be one of: "
180 >                "MIX, MAX, or TRADITIONAL\n",
181 >                cutPolicy.c_str());
182 >        painCave.isFatal = 1;
183 >        painCave.severity = OPENMD_ERROR;
184 >        simError();
185 >      } else {
186 >        cutoffPolicy_ = i->second;
187 >      }
188 >    } else {
189 >      sprintf(painCave.errMsg,
190 >              "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
191 >              "\tOpenMD will use TRADITIONAL.\n");
192 >      painCave.isFatal = 0;
193 >      painCave.severity = OPENMD_INFO;
194 >      simError();
195 >      cutoffPolicy_ = TRADITIONAL;        
196 >    }
197 >  }
198 >
199 >  /**
200 >   * setupSwitching
201 >   *
202 >   * Sets the values of switchingRadius and
203 >   *  If the switchingRadius was explicitly set, use that value (but check it)
204 >   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
205 >   */
206 >  void ForceManager::setupSwitching() {
207 >    Globals* simParams_ = info_->getSimParams();
208      
209 <    if (!info_->isFortranInitialized()) {
209 >    if (simParams_->haveSwitchingRadius()) {
210 >      rSwitch_ = simParams_->getSwitchingRadius();
211 >      if (rSwitch_ > rCut_) {        
212 >        sprintf(painCave.errMsg,
213 >                "ForceManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n",
214 >                rSwitch_, rCut_);
215 >        painCave.isFatal = 1;
216 >        painCave.severity = OPENMD_ERROR;
217 >        simError();
218 >      }
219 >    } else {      
220 >      rSwitch_ = 0.85 * rCut_;
221 >      sprintf(painCave.errMsg,
222 >              "ForceManager::setupSwitching: No value was set for the switchingRadius.\n"
223 >              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
224 >              "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
225 >      painCave.isFatal = 0;
226 >      painCave.severity = OPENMD_WARNING;
227 >      simError();
228 >    }          
229 >    
230 >    if (simParams_->haveSwitchingFunctionType()) {
231 >      string funcType = simParams_->getSwitchingFunctionType();
232 >      toUpper(funcType);
233 >      if (funcType == "CUBIC") {
234 >        sft_ = cubic;
235 >      } else {
236 >        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
237 >          sft_ = fifth_order_poly;
238 >        } else {
239 >          // throw error        
240 >          sprintf( painCave.errMsg,
241 >                   "ForceManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n"
242 >                   "\tswitchingFunctionType must be one of: "
243 >                   "\"cubic\" or \"fifth_order_polynomial\".",
244 >                   funcType.c_str() );
245 >          painCave.isFatal = 1;
246 >          painCave.severity = OPENMD_ERROR;
247 >          simError();
248 >        }          
249 >      }
250 >    }
251 >    switcher_->setSwitchType(sft_);
252 >    switcher_->setSwitch(rSwitch_, rCut_);
253 >  }
254 >  
255 >  void ForceManager::initialize() {
256 >
257 >    if (!info_->isTopologyDone()) {
258        info_->update();
259 <      nbiMan_->setSimInfo(info_);
260 <      nbiMan_->initialize();
261 <      swfun_ = nbiMan_->getSwitchingFunction();
262 <      decomp_->distributeInitialData();
263 <      info_->setupFortran();
259 >      interactionMan_->setSimInfo(info_);
260 >      interactionMan_->initialize();
261 >
262 >      // We want to delay the cutoffs until after the interaction
263 >      // manager has set up the atom-atom interactions so that we can
264 >      // query them for suggested cutoff values
265 >
266 >      setupCutoffs();
267 >      setupSwitching();
268 >
269 >      info_->prepareTopology();      
270      }
271 +
272 +    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
273      
274 <    preCalculation();  
275 <    calcShortRangeInteraction();
276 <    calcLongRangeInteraction();
277 <    postCalculation();
274 >    // Force fields can set options on how to scale van der Waals and electrostatic
275 >    // interactions for atoms connected via bonds, bends and torsions
276 >    // in this case the topological distance between atoms is:
277 >    // 0 = topologically unconnected
278 >    // 1 = bonded together
279 >    // 2 = connected via a bend
280 >    // 3 = connected via a torsion
281 >    
282 >    vdwScale_.reserve(4);
283 >    fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
284 >
285 >    electrostaticScale_.reserve(4);
286 >    fill(electrostaticScale_.begin(), electrostaticScale_.end(), 0.0);
287 >
288 >    vdwScale_[0] = 1.0;
289 >    vdwScale_[1] = fopts.getvdw12scale();
290 >    vdwScale_[2] = fopts.getvdw13scale();
291 >    vdwScale_[3] = fopts.getvdw14scale();
292 >    
293 >    electrostaticScale_[0] = 1.0;
294 >    electrostaticScale_[1] = fopts.getelectrostatic12scale();
295 >    electrostaticScale_[2] = fopts.getelectrostatic13scale();
296 >    electrostaticScale_[3] = fopts.getelectrostatic14scale();    
297 >    
298 >    fDecomp_->distributeInitialData();
299 >
300 >    initialized_ = true;
301 >
302 >  }
303 >
304 >  void ForceManager::calcForces() {
305      
306 +    if (!initialized_) initialize();
307 +
308 +    preCalculation();  
309 +    shortRangeInteractions();
310 +    longRangeInteractions();
311 +    postCalculation();    
312    }
313    
314    void ForceManager::preCalculation() {
# Line 128 | Line 349 | namespace OpenMD {
349      
350    }
351    
352 <  void ForceManager::calcShortRangeInteraction() {
352 >  void ForceManager::shortRangeInteractions() {
353      Molecule* mol;
354      RigidBody* rb;
355      Bond* bond;
# Line 245 | Line 466 | namespace OpenMD {
466      curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
467    }
468    
469 <  void ForceManager::calcLongRangeInteraction() {
469 >  void ForceManager::longRangeInteractions() {
470  
471      // some of this initial stuff will go away:
472      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 267 | Line 488 | namespace OpenMD {
488        rc = pos;
489      }
490      
270 –    //initialize data before passing to fortran
271 –    RealType longRangePotential[LR_POT_TYPES];
272 –    RealType lrPot = 0.0;
273 –    int isError = 0;
274 –
275 –    for (int i=0; i<LR_POT_TYPES;i++){
276 –      longRangePotential[i]=0.0; //Initialize array
277 –    }
278 –
491      // new stuff starts here:
492 <
493 <    decomp_->distributeData();
492 >    fDecomp_->zeroWorkArrays();
493 >    fDecomp_->distributeData();
494  
495 <    int cg1, cg2;
496 <    Vector3d d_grp;
495 >    int cg1, cg2, atom1, atom2;
496 >    Vector3d d_grp, dag;
497      RealType rgrpsq, rgrp;
498      RealType vij;
499      Vector3d fij, fg;
500 <    pair<int, int> gtypes;
500 >    tuple3<RealType, RealType, RealType> cuts;
501      RealType rCutSq;
502      bool in_switching_region;
503      RealType sw, dswdr, swderiv;
504 <    vector<int> atomListI;
293 <    vector<int> atomListJ;
504 >    vector<int> atomListColumn, atomListRow, atomListLocal;
505      InteractionData idat;
506 +    SelfData sdat;
507 +    RealType mf;
508 +    potVec pot(0.0);
509 +    potVec longRangePotential(0.0);
510 +    RealType lrPot;
511  
512      int loopStart, loopEnd;
513  
514      loopEnd = PAIR_LOOP;
515 <    if (info_->requiresPrepair_) {
515 >    if (info_->requiresPrepair() ) {
516        loopStart = PREPAIR_LOOP;
517      } else {
518        loopStart = PAIR_LOOP;
# Line 305 | Line 521 | namespace OpenMD {
521      for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
522        
523        if (iLoop == loopStart) {
524 <        bool update_nlist = decomp_->checkNeighborList();
524 >        bool update_nlist = fDecomp_->checkNeighborList();
525          if (update_nlist)
526 <          neighborList = decomp_->buildNeighborList();
526 >          neighborList = fDecomp_->buildNeighborList();
527        }
528  
529        for (vector<pair<int, int> >::iterator it = neighborList.begin();
# Line 315 | Line 531 | namespace OpenMD {
531          
532          cg1 = (*it).first;
533          cg2 = (*it).second;
534 +        
535 +        cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
536  
537 <        gtypes = decomp_->getGroupTypes(cg1, cg2);
320 <        d_grp  = decomp_->getIntergroupVector(cg1, cg2);
537 >        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
538          curSnapshot->wrapVector(d_grp);        
539          rgrpsq = d_grp.lengthSquare();
323 –        rCutSq = groupCutoffMap(gtypes).first;
540  
541 +        rCutSq = cuts.second;
542 +
543          if (rgrpsq < rCutSq) {
544 <          idat.rcut = groupCutoffMap(gtypes).second;
544 >          *(idat.rcut) = cuts.first;
545            if (iLoop == PAIR_LOOP) {
546 <            vij = 0.0;
546 >            vij *= 0.0;
547              fij = V3Zero;
548            }
549            
550 <          in_switching_region = swfun_->getSwitch(rgrpsq, idat.sw, idat.dswdr, rgrp);    
551 <          
552 <          atomListI = decomp_->getAtomsInGroupI(cg1);
553 <          atomListJ = decomp_->getAtomsInGroupJ(cg2);
550 >          in_switching_region = switcher_->getSwitch(rgrpsq, *(idat.sw), dswdr,
551 >                                                     rgrp);
552 >              
553 >          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
554 >          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
555  
556 <          for (vector<int>::iterator ia = atomListI.begin();
557 <               ia != atomListI.end(); ++ia) {            
556 >          for (vector<int>::iterator ia = atomListRow.begin();
557 >               ia != atomListRow.end(); ++ia) {            
558              atom1 = (*ia);
559              
560 <            for (vector<int>::iterator jb = atomListJ.begin();
561 <                 jb != atomListJ.end(); ++jb) {              
560 >            for (vector<int>::iterator jb = atomListColumn.begin();
561 >                 jb != atomListColumn.end(); ++jb) {              
562                atom2 = (*jb);
563                
564 <              if (!decomp_->skipAtomPair(atom1, atom2)) {
564 >              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
565                  
566 <                if (atomListI.size() == 1 && atomListJ.size() == 1) {
567 <                  idat.d = d_grp;
568 <                  idat.r2 = rgrpsq;
566 >                pot *= 0.0;
567 >
568 >                idat = fDecomp_->fillInteractionData(atom1, atom2);
569 >                *(idat.pot) = pot;
570 >
571 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
572 >                  *(idat.d) = d_grp;
573 >                  *(idat.r2) = rgrpsq;
574                  } else {
575 <                  idat.d = decomp_->getInteratomicVector(atom1, atom2);
576 <                  curSnapshot->wrapVector(idat.d);
577 <                  idat.r2 = idat.d.lengthSquare();
575 >                  *(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2);
576 >                  curSnapshot->wrapVector( *(idat.d) );
577 >                  *(idat.r2) = idat.d->lengthSquare();
578                  }
579                  
580 <                idat.r = sqrt(idat.r2);
581 <                decomp_->fillInteractionData(atom1, atom2, idat);
358 <                
580 >                *(idat.rij) = sqrt( *(idat.r2) );
581 >              
582                  if (iLoop == PREPAIR_LOOP) {
583                    interactionMan_->doPrePair(idat);
584                  } else {
585                    interactionMan_->doPair(idat);
586 <                  vij += idat.vpair;
587 <                  fij += idat.f1;
588 <                  tau -= outProduct(idat.d, idat.f);
586 >                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
587 >                  vij += *(idat.vpair);
588 >                  fij += *(idat.f1);
589 >                  tau -= outProduct( *(idat.d), *(idat.f1));
590                  }
591                }
592              }
# Line 375 | Line 599 | namespace OpenMD {
599  
600                fij += fg;
601  
602 <              if (atomListI.size() == 1 && atomListJ.size() == 1) {
603 <                tau -= outProduct(idat.d, fg);
602 >              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
603 >                tau -= outProduct( *(idat.d), fg);
604                }
605            
606 <              for (vector<int>::iterator ia = atomListI.begin();
607 <                   ia != atomListI.end(); ++ia) {            
606 >              for (vector<int>::iterator ia = atomListRow.begin();
607 >                   ia != atomListRow.end(); ++ia) {            
608                  atom1 = (*ia);                
609 <                mf = decomp_->getMfactI(atom1);
609 >                mf = fDecomp_->getMassFactorRow(atom1);
610                  // fg is the force on atom ia due to cutoff group's
611                  // presence in switching region
612                  fg = swderiv * d_grp * mf;
613 <                decomp_->addForceToAtomI(atom1, fg);
613 >                fDecomp_->addForceToAtomRow(atom1, fg);
614  
615 <                if (atomListI.size() > 1) {
616 <                  if (info_->usesAtomicVirial_) {
615 >                if (atomListRow.size() > 1) {
616 >                  if (info_->usesAtomicVirial()) {
617                      // find the distance between the atom
618                      // and the center of the cutoff group:
619 <                    dag = decomp_->getAtomToGroupVectorI(atom1, cg1);
619 >                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
620                      tau -= outProduct(dag, fg);
621                    }
622                  }
623                }
624 <              for (vector<int>::iterator jb = atomListJ.begin();
625 <                   jb != atomListJ.end(); ++jb) {              
624 >              for (vector<int>::iterator jb = atomListColumn.begin();
625 >                   jb != atomListColumn.end(); ++jb) {              
626                  atom2 = (*jb);
627 <                mf = decomp_->getMfactJ(atom2);
627 >                mf = fDecomp_->getMassFactorColumn(atom2);
628                  // fg is the force on atom jb due to cutoff group's
629                  // presence in switching region
630                  fg = -swderiv * d_grp * mf;
631 <                decomp_->addForceToAtomJ(atom2, fg);
631 >                fDecomp_->addForceToAtomColumn(atom2, fg);
632  
633 <                if (atomListJ.size() > 1) {
634 <                  if (info_->usesAtomicVirial_) {
633 >                if (atomListColumn.size() > 1) {
634 >                  if (info_->usesAtomicVirial()) {
635                      // find the distance between the atom
636                      // and the center of the cutoff group:
637 <                    dag = decomp_->getAtomToGroupVectorJ(atom2, cg2);
637 >                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
638                      tau -= outProduct(dag, fg);
639                    }
640                  }
# Line 424 | Line 648 | namespace OpenMD {
648        }
649  
650        if (iLoop == PREPAIR_LOOP) {
651 <        if (info_->requiresPrepair_) {            
652 <          decomp_->collectIntermediateData();
653 <          atomList = decomp_->getAtomList();
654 <          for (vector<int>::iterator ia = atomList.begin();
655 <               ia != atomList.end(); ++ia) {              
432 <            atom1 = (*ia);            
433 <            decomp_->populateSelfData(atom1, SelfData sdat);
651 >        if (info_->requiresPrepair()) {            
652 >          fDecomp_->collectIntermediateData();
653 >
654 >          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
655 >            sdat = fDecomp_->fillSelfData(atom1);
656              interactionMan_->doPreForce(sdat);
657            }
658 <          decomp_->distributeIntermediateData();        
658 >
659 >          fDecomp_->distributeIntermediateData();        
660          }
661        }
662  
663      }
664      
665 <    decomp_->collectData();
665 >    fDecomp_->collectData();
666      
667 <    if (info_->requiresSkipCorrection_ || info_->requiresSelfCorrection_) {
668 <      atomList = decomp_->getAtomList();
669 <      for (vector<int>::iterator ia = atomList.begin();
447 <           ia != atomList.end(); ++ia) {              
448 <        atom1 = (*ia);    
667 >    if ( info_->requiresSkipCorrection() ) {
668 >      
669 >      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
670  
671 <        if (info_->requiresSkipCorrection_) {
672 <          vector<int> skipList = decomp_->getSkipsForAtom(atom1);
673 <          for (vector<int>::iterator jb = skipList.begin();
674 <               jb != skipList.end(); ++jb) {              
675 <            atom2 = (*jb);
676 <            decomp_->populateSkipData(atom1, atom2, InteractionData idat);
677 <            interactionMan_->doSkipCorrection(idat);
678 <          }
671 >        vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
672 >        
673 >        for (vector<int>::iterator jb = skipList.begin();
674 >             jb != skipList.end(); ++jb) {        
675 >    
676 >          atom2 = (*jb);
677 >          idat = fDecomp_->fillSkipData(atom1, atom2);
678 >          interactionMan_->doSkipCorrection(idat);
679 >
680          }
459 –          
460 –        if (info_->requiresSelfCorrection_) {
461 –          decomp_->populateSelfData(atom1, SelfData sdat);
462 –          interactionMan_->doSelfCorrection(sdat);
681        }
464 –      
465 –      
682      }
683 +    
684 +    if (info_->requiresSelfCorrection()) {
685  
686 <    for (int i=0; i<LR_POT_TYPES;i++){
687 <      lrPot += longRangePotential[i]; //Quick hack
686 >      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
687 >        sdat = fDecomp_->fillSelfData(atom1);
688 >        interactionMan_->doSelfCorrection(sdat);
689 >      }
690 >
691      }
692 <        
692 >
693 >    longRangePotential = fDecomp_->getLongRangePotential();
694 >    lrPot = longRangePotential.sum();
695 >
696      //store the tau and long range potential    
697      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
698 <    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
699 <    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
698 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
699 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
700    }
701  
702    

Diff Legend

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