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 1470 by gezelter, Mon Jul 19 18:48:23 2010 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 "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 <                                               NBforcesInitialized_(false) {
67 <    lj_ = LJ::Instance();
66 <    lj_->setForceField(info_->getForceField());
65 >  ForceManager::ForceManager(SimInfo * info) : info_(info) {
66 >    forceField_ = info_->getForceField();
67 >    fDecomp_ = new ForceMatrixDecomposition(info_);
68    }
69 <
70 <  void ForceManager::calcForces() {
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 <    if (!info_->isFortranInitialized()) {
94 <      info_->update();
73 <    }
93 >    Globals* simParams_ = info_->getSimParams();
94 >    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
95      
96 <    preCalculation();
97 <    
98 <    calcShortRangeInteraction();
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 <    calcLongRangeInteraction();
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 >    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 <    postCalculation();
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 (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 +      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 +    // 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() {
315      SimInfo::MoleculeIterator mi;
316      Molecule* mol;
# Line 89 | Line 318 | namespace OpenMD {
318      Atom* atom;
319      Molecule::RigidBodyIterator rbIter;
320      RigidBody* rb;
321 +    Molecule::CutoffGroupIterator ci;
322 +    CutoffGroup* cg;
323      
324      // forces are zeroed here, before any are accumulated.
94 –    // NOTE: do not rezero the forces in Fortran.
325      
326      for (mol = info_->beginMolecule(mi); mol != NULL;
327           mol = info_->nextMolecule(mi)) {
# Line 104 | Line 334 | namespace OpenMD {
334             rb = mol->nextRigidBody(rbIter)) {
335          rb->zeroForcesAndTorques();
336        }        
337 <          
337 >
338 >      if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
339 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
340 >            cg = mol->nextCutoffGroup(ci)) {
341 >          //calculate the center of mass of cutoff group
342 >          cg->updateCOM();
343 >        }
344 >      }      
345      }
346 <    
346 >  
347      // Zero out the stress tensor
348      tau *= 0.0;
349      
350    }
351    
352 <  void ForceManager::calcShortRangeInteraction() {
352 >  void ForceManager::shortRangeInteractions() {
353      Molecule* mol;
354      RigidBody* rb;
355      Bond* bond;
# Line 154 | Line 391 | namespace OpenMD {
391          RealType currBendPot = bend->getPotential();          
392          
393          bendPotential += bend->getPotential();
394 <        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
394 >        map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
395          if (i == bendDataSets.end()) {
396            BendDataSet dataSet;
397            dataSet.prev.angle = dataSet.curr.angle = angle;
398            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
399            dataSet.deltaV = 0.0;
400 <          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
400 >          bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
401          }else {
402            i->second.prev.angle = i->second.curr.angle;
403            i->second.prev.potential = i->second.curr.potential;
# Line 177 | Line 414 | namespace OpenMD {
414          torsion->calcForce(angle);
415          RealType currTorsionPot = torsion->getPotential();
416          torsionPotential += torsion->getPotential();
417 <        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
417 >        map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
418          if (i == torsionDataSets.end()) {
419            TorsionDataSet dataSet;
420            dataSet.prev.angle = dataSet.curr.angle = angle;
421            dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
422            dataSet.deltaV = 0.0;
423 <          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
423 >          torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
424          }else {
425            i->second.prev.angle = i->second.curr.angle;
426            i->second.prev.potential = i->second.curr.potential;
# Line 193 | Line 430 | namespace OpenMD {
430                                     i->second.prev.potential);
431          }      
432        }      
433 <
433 >      
434        for (inversion = mol->beginInversion(inversionIter);
435             inversion != NULL;
436             inversion = mol->nextInversion(inversionIter)) {
# Line 201 | Line 438 | namespace OpenMD {
438          inversion->calcForce(angle);
439          RealType currInversionPot = inversion->getPotential();
440          inversionPotential += inversion->getPotential();
441 <        std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
441 >        map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
442          if (i == inversionDataSets.end()) {
443            InversionDataSet dataSet;
444            dataSet.prev.angle = dataSet.curr.angle = angle;
445            dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
446            dataSet.deltaV = 0.0;
447 <          inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
447 >          inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
448          }else {
449            i->second.prev.angle = i->second.curr.angle;
450            i->second.prev.potential = i->second.curr.potential;
# Line 226 | Line 463 | namespace OpenMD {
463      curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
464      curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
465      curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
466 <    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
230 <    
466 >    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
467    }
468    
469 <  void ForceManager::calcLongRangeInteraction() {
234 <    Snapshot* curSnapshot;
235 <    DataStorage* config;
236 <    RealType* frc;
237 <    RealType* pos;
238 <    RealType* trq;
239 <    RealType* A;
240 <    RealType* electroFrame;
241 <    RealType* rc;
242 <    RealType* particlePot;
243 <    
244 <    //get current snapshot from SimInfo
245 <    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
246 <    
247 <    //get array pointers
248 <    config = &(curSnapshot->atomData);
249 <    frc = config->getArrayPointer(DataStorage::dslForce);
250 <    pos = config->getArrayPointer(DataStorage::dslPosition);
251 <    trq = config->getArrayPointer(DataStorage::dslTorque);
252 <    A   = config->getArrayPointer(DataStorage::dslAmat);
253 <    electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
254 <    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
469 >  void ForceManager::longRangeInteractions() {
470  
471 <    //calculate the center of mass of cutoff group
472 <    SimInfo::MoleculeIterator mi;
473 <    Molecule* mol;
474 <    Molecule::CutoffGroupIterator ci;
475 <    CutoffGroup* cg;
476 <    Vector3d com;
477 <    std::vector<Vector3d> rcGroup;
478 <    
479 <    if(info_->getNCutoffGroups() > 0){
480 <      
481 <      for (mol = info_->beginMolecule(mi); mol != NULL;
482 <           mol = info_->nextMolecule(mi)) {
483 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
484 <            cg = mol->nextCutoffGroup(ci)) {
270 <          cg->getCOM(com);
271 <          rcGroup.push_back(com);
272 <        }
273 <      }// end for (mol)
274 <      
275 <      rc = rcGroup[0].getArrayPointer();
471 >    // some of this initial stuff will go away:
472 >    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
473 >    DataStorage* config = &(curSnapshot->atomData);
474 >    DataStorage* cgConfig = &(curSnapshot->cgData);
475 >    RealType* frc = config->getArrayPointer(DataStorage::dslForce);
476 >    RealType* pos = config->getArrayPointer(DataStorage::dslPosition);
477 >    RealType* trq = config->getArrayPointer(DataStorage::dslTorque);
478 >    RealType* A = config->getArrayPointer(DataStorage::dslAmat);
479 >    RealType* electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
480 >    RealType* particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
481 >    RealType* rc;    
482 >
483 >    if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
484 >      rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
485      } else {
486        // center of mass of the group is the same as position of the atom  
487        // if cutoff group does not exist
488        rc = pos;
489      }
490      
491 <    //initialize data before passing to fortran
492 <    RealType longRangePotential[LR_POT_TYPES];
493 <    RealType lrPot = 0.0;
494 <    Vector3d totalDipole;
495 <    int isError = 0;
491 >    // new stuff starts here:
492 >    fDecomp_->zeroWorkArrays();
493 >    fDecomp_->distributeData();
494 >
495 >    int cg1, cg2, atom1, atom2;
496 >    Vector3d d_grp, dag;
497 >    RealType rgrpsq, rgrp;
498 >    RealType vij;
499 >    Vector3d fij, fg;
500 >    tuple3<RealType, RealType, RealType> cuts;
501 >    RealType rCutSq;
502 >    bool in_switching_region;
503 >    RealType sw, dswdr, swderiv;
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 <    for (int i=0; i<LR_POT_TYPES;i++){
513 <      longRangePotential[i]=0.0; //Initialize array
512 >    int loopStart, loopEnd;
513 >
514 >    loopEnd = PAIR_LOOP;
515 >    if (info_->requiresPrepair() ) {
516 >      loopStart = PREPAIR_LOOP;
517 >    } else {
518 >      loopStart = PAIR_LOOP;
519      }
520 <    
521 <    doForceLoop(pos,
522 <                rc,
523 <                A,
524 <                electroFrame,
525 <                frc,
526 <                trq,
527 <                tau.getArrayPointer(),
528 <                longRangePotential,
529 <                particlePot,
530 <                &isError );
531 <    
532 <    if( isError ){
533 <      sprintf( painCave.errMsg,
534 <               "Error returned from the fortran force calculation.\n" );
535 <      painCave.isFatal = 1;
536 <      simError();
520 >
521 >    for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
522 >      
523 >      if (iLoop == loopStart) {
524 >        bool update_nlist = fDecomp_->checkNeighborList();
525 >        if (update_nlist)
526 >          neighborList = fDecomp_->buildNeighborList();
527 >      }
528 >
529 >      for (vector<pair<int, int> >::iterator it = neighborList.begin();
530 >             it != neighborList.end(); ++it) {
531 >        
532 >        cg1 = (*it).first;
533 >        cg2 = (*it).second;
534 >        
535 >        cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
536 >
537 >        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
538 >        curSnapshot->wrapVector(d_grp);        
539 >        rgrpsq = d_grp.lengthSquare();
540 >
541 >        rCutSq = cuts.second;
542 >
543 >        if (rgrpsq < rCutSq) {
544 >          *(idat.rcut) = cuts.first;
545 >          if (iLoop == PAIR_LOOP) {
546 >            vij *= 0.0;
547 >            fij = V3Zero;
548 >          }
549 >          
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 = atomListRow.begin();
557 >               ia != atomListRow.end(); ++ia) {            
558 >            atom1 = (*ia);
559 >            
560 >            for (vector<int>::iterator jb = atomListColumn.begin();
561 >                 jb != atomListColumn.end(); ++jb) {              
562 >              atom2 = (*jb);
563 >              
564 >              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
565 >                
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) = fDecomp_->getInteratomicVector(atom1, atom2);
576 >                  curSnapshot->wrapVector( *(idat.d) );
577 >                  *(idat.r2) = idat.d->lengthSquare();
578 >                }
579 >                
580 >                *(idat.rij) = sqrt( *(idat.r2) );
581 >              
582 >                if (iLoop == PREPAIR_LOOP) {
583 >                  interactionMan_->doPrePair(idat);
584 >                } else {
585 >                  interactionMan_->doPair(idat);
586 >                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
587 >                  vij += *(idat.vpair);
588 >                  fij += *(idat.f1);
589 >                  tau -= outProduct( *(idat.d), *(idat.f1));
590 >                }
591 >              }
592 >            }
593 >          }
594 >
595 >          if (iLoop == PAIR_LOOP) {
596 >            if (in_switching_region) {
597 >              swderiv = vij * dswdr / rgrp;
598 >              fg = swderiv * d_grp;
599 >
600 >              fij += fg;
601 >
602 >              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
603 >                tau -= outProduct( *(idat.d), fg);
604 >              }
605 >          
606 >              for (vector<int>::iterator ia = atomListRow.begin();
607 >                   ia != atomListRow.end(); ++ia) {            
608 >                atom1 = (*ia);                
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 >                fDecomp_->addForceToAtomRow(atom1, fg);
614 >
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 = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
620 >                    tau -= outProduct(dag, fg);
621 >                  }
622 >                }
623 >              }
624 >              for (vector<int>::iterator jb = atomListColumn.begin();
625 >                   jb != atomListColumn.end(); ++jb) {              
626 >                atom2 = (*jb);
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 >                fDecomp_->addForceToAtomColumn(atom2, fg);
632 >
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 = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
638 >                    tau -= outProduct(dag, fg);
639 >                  }
640 >                }
641 >              }
642 >            }
643 >            //if (!SIM_uses_AtomicVirial) {
644 >            //  tau -= outProduct(d_grp, fij);
645 >            //}
646 >          }
647 >        }
648 >      }
649 >
650 >      if (iLoop == PREPAIR_LOOP) {
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 >
659 >          fDecomp_->distributeIntermediateData();        
660 >        }
661 >      }
662 >
663      }
309 –    for (int i=0; i<LR_POT_TYPES;i++){
310 –      lrPot += longRangePotential[i]; //Quick hack
311 –    }
664      
665 <    // grab the simulation box dipole moment if specified
666 <    if (info_->getCalcBoxDipole()){
667 <      getAccumulatedBoxDipole(totalDipole.getArrayPointer());
665 >    fDecomp_->collectData();
666 >    
667 >    if ( info_->requiresSkipCorrection() ) {
668        
669 <      curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
670 <      curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
671 <      curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
669 >      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
670 >
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 >        }
681 >      }
682      }
683      
684 +    if (info_->requiresSelfCorrection()) {
685 +
686 +      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
687 +        sdat = fDecomp_->fillSelfData(atom1);
688 +        interactionMan_->doSelfCorrection(sdat);
689 +      }
690 +
691 +    }
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