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

Diff Legend

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