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 1503 by gezelter, Sat Oct 2 19:54:41 2010 UTC vs.
Revision 1584 by gezelter, Fri Jun 17 20:16:35 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 + #include <cstdio>
63 + #include <iostream>
64 + #include <iomanip>
65 +
66 + using namespace std;
67   namespace OpenMD {
68    
69 <  ForceManager::ForceManager(SimInfo * info) : info_(info),
70 <                                               NBforcesInitialized_(false) {
69 >  ForceManager::ForceManager(SimInfo * info) : info_(info) {
70 >    forceField_ = info_->getForceField();
71 >    interactionMan_ = new InteractionManager();
72 >    fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_);
73    }
74 <
75 <  void ForceManager::calcForces() {
74 >
75 >  /**
76 >   * setupCutoffs
77 >   *
78 >   * Sets the values of cutoffRadius, cutoffMethod, and cutoffPolicy
79 >   *
80 >   * cutoffRadius : realType
81 >   *  If the cutoffRadius was explicitly set, use that value.
82 >   *  If the cutoffRadius was not explicitly set:
83 >   *      Are there electrostatic atoms?  Use 12.0 Angstroms.
84 >   *      No electrostatic atoms?  Poll the atom types present in the
85 >   *      simulation for suggested cutoff values (e.g. 2.5 * sigma).
86 >   *      Use the maximum suggested value that was found.
87 >   *
88 >   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
89 >   *      If cutoffMethod was explicitly set, use that choice.
90 >   *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
91 >   *
92 >   * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
93 >   *      If cutoffPolicy was explicitly set, use that choice.
94 >   *      If cutoffPolicy was not explicitly set, use TRADITIONAL
95 >   */
96 >  void ForceManager::setupCutoffs() {
97      
98 <    if (!info_->isFortranInitialized()) {
99 <      info_->update();
71 <    }
98 >    Globals* simParams_ = info_->getSimParams();
99 >    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
100      
101 <    preCalculation();
102 <    
103 <    calcShortRangeInteraction();
101 >    if (simParams_->haveCutoffRadius()) {
102 >      rCut_ = simParams_->getCutoffRadius();
103 >    } else {      
104 >      if (info_->usesElectrostaticAtoms()) {
105 >        sprintf(painCave.errMsg,
106 >                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
107 >                "\tOpenMD will use a default value of 12.0 angstroms"
108 >                "\tfor the cutoffRadius.\n");
109 >        painCave.isFatal = 0;
110 >        painCave.severity = OPENMD_INFO;
111 >        simError();
112 >        rCut_ = 12.0;
113 >      } else {
114 >        RealType thisCut;
115 >        set<AtomType*>::iterator i;
116 >        set<AtomType*> atomTypes;
117 >        atomTypes = info_->getSimulatedAtomTypes();        
118 >        for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
119 >          thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
120 >          rCut_ = max(thisCut, rCut_);
121 >        }
122 >        sprintf(painCave.errMsg,
123 >                "ForceManager::setupCutoffs: No value was set for the cutoffRadius.\n"
124 >                "\tOpenMD will use %lf angstroms.\n",
125 >                rCut_);
126 >        painCave.isFatal = 0;
127 >        painCave.severity = OPENMD_INFO;
128 >        simError();
129 >      }
130 >    }
131  
132 <    calcLongRangeInteraction();
132 >    fDecomp_->setUserCutoff(rCut_);
133 >    interactionMan_->setCutoffRadius(rCut_);
134  
135 <    postCalculation();
135 >    map<string, CutoffMethod> stringToCutoffMethod;
136 >    stringToCutoffMethod["HARD"] = HARD;
137 >    stringToCutoffMethod["SWITCHED"] = SWITCHED;
138 >    stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;    
139 >    stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
140 >  
141 >    if (simParams_->haveCutoffMethod()) {
142 >      string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
143 >      map<string, CutoffMethod>::iterator i;
144 >      i = stringToCutoffMethod.find(cutMeth);
145 >      if (i == stringToCutoffMethod.end()) {
146 >        sprintf(painCave.errMsg,
147 >                "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
148 >                "\tShould be one of: "
149 >                "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
150 >                cutMeth.c_str());
151 >        painCave.isFatal = 1;
152 >        painCave.severity = OPENMD_ERROR;
153 >        simError();
154 >      } else {
155 >        cutoffMethod_ = i->second;
156 >      }
157 >    } else {
158 >      sprintf(painCave.errMsg,
159 >              "ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n"
160 >              "\tOpenMD will use SHIFTED_FORCE.\n");
161 >      painCave.isFatal = 0;
162 >      painCave.severity = OPENMD_INFO;
163 >      simError();
164 >      cutoffMethod_ = SHIFTED_FORCE;        
165 >    }
166 >
167 >    map<string, CutoffPolicy> stringToCutoffPolicy;
168 >    stringToCutoffPolicy["MIX"] = MIX;
169 >    stringToCutoffPolicy["MAX"] = MAX;
170 >    stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;    
171 >
172 >    std::string cutPolicy;
173 >    if (forceFieldOptions_.haveCutoffPolicy()){
174 >      cutPolicy = forceFieldOptions_.getCutoffPolicy();
175 >    }else if (simParams_->haveCutoffPolicy()) {
176 >      cutPolicy = simParams_->getCutoffPolicy();
177 >    }
178 >
179 >    if (!cutPolicy.empty()){
180 >      toUpper(cutPolicy);
181 >      map<string, CutoffPolicy>::iterator i;
182 >      i = stringToCutoffPolicy.find(cutPolicy);
183 >
184 >      if (i == stringToCutoffPolicy.end()) {
185 >        sprintf(painCave.errMsg,
186 >                "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
187 >                "\tShould be one of: "
188 >                "MIX, MAX, or TRADITIONAL\n",
189 >                cutPolicy.c_str());
190 >        painCave.isFatal = 1;
191 >        painCave.severity = OPENMD_ERROR;
192 >        simError();
193 >      } else {
194 >        cutoffPolicy_ = i->second;
195 >      }
196 >    } else {
197 >      sprintf(painCave.errMsg,
198 >              "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
199 >              "\tOpenMD will use TRADITIONAL.\n");
200 >      painCave.isFatal = 0;
201 >      painCave.severity = OPENMD_INFO;
202 >      simError();
203 >      cutoffPolicy_ = TRADITIONAL;        
204 >    }
205 >    fDecomp_->setCutoffPolicy(cutoffPolicy_);
206 >  }
207 >
208 >  /**
209 >   * setupSwitching
210 >   *
211 >   * Sets the values of switchingRadius and
212 >   *  If the switchingRadius was explicitly set, use that value (but check it)
213 >   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
214 >   */
215 >  void ForceManager::setupSwitching() {
216 >    Globals* simParams_ = info_->getSimParams();
217 >
218 >    // create the switching function object:
219 >    switcher_ = new SwitchingFunction();
220      
221 +    if (simParams_->haveSwitchingRadius()) {
222 +      rSwitch_ = simParams_->getSwitchingRadius();
223 +      if (rSwitch_ > rCut_) {        
224 +        sprintf(painCave.errMsg,
225 +                "ForceManager::setupSwitching: switchingRadius (%f) is larger "
226 +                "than the cutoffRadius(%f)\n", rSwitch_, rCut_);
227 +        painCave.isFatal = 1;
228 +        painCave.severity = OPENMD_ERROR;
229 +        simError();
230 +      }
231 +    } else {      
232 +      rSwitch_ = 0.85 * rCut_;
233 +      sprintf(painCave.errMsg,
234 +              "ForceManager::setupSwitching: No value was set for the switchingRadius.\n"
235 +              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
236 +              "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
237 +      painCave.isFatal = 0;
238 +      painCave.severity = OPENMD_WARNING;
239 +      simError();
240 +    }          
241 +    
242 +    // Default to cubic switching function.
243 +    sft_ = cubic;
244 +    if (simParams_->haveSwitchingFunctionType()) {
245 +      string funcType = simParams_->getSwitchingFunctionType();
246 +      toUpper(funcType);
247 +      if (funcType == "CUBIC") {
248 +        sft_ = cubic;
249 +      } else {
250 +        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
251 +          sft_ = fifth_order_poly;
252 +        } else {
253 +          // throw error        
254 +          sprintf( painCave.errMsg,
255 +                   "ForceManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n"
256 +                   "\tswitchingFunctionType must be one of: "
257 +                   "\"cubic\" or \"fifth_order_polynomial\".",
258 +                   funcType.c_str() );
259 +          painCave.isFatal = 1;
260 +          painCave.severity = OPENMD_ERROR;
261 +          simError();
262 +        }          
263 +      }
264 +    }
265 +    switcher_->setSwitchType(sft_);
266 +    switcher_->setSwitch(rSwitch_, rCut_);
267 +    interactionMan_->setSwitchingRadius(rSwitch_);
268    }
269    
270 +  void ForceManager::initialize() {
271 +
272 +    if (!info_->isTopologyDone()) {
273 +      info_->update();
274 +      interactionMan_->setSimInfo(info_);
275 +      interactionMan_->initialize();
276 +
277 +      // We want to delay the cutoffs until after the interaction
278 +      // manager has set up the atom-atom interactions so that we can
279 +      // query them for suggested cutoff values
280 +
281 +      setupCutoffs();
282 +      setupSwitching();
283 +
284 +      info_->prepareTopology();      
285 +    }
286 +
287 +    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
288 +    
289 +    // Force fields can set options on how to scale van der Waals and electrostatic
290 +    // interactions for atoms connected via bonds, bends and torsions
291 +    // in this case the topological distance between atoms is:
292 +    // 0 = topologically unconnected
293 +    // 1 = bonded together
294 +    // 2 = connected via a bend
295 +    // 3 = connected via a torsion
296 +    
297 +    vdwScale_.reserve(4);
298 +    fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
299 +
300 +    electrostaticScale_.reserve(4);
301 +    fill(electrostaticScale_.begin(), electrostaticScale_.end(), 0.0);
302 +
303 +    vdwScale_[0] = 1.0;
304 +    vdwScale_[1] = fopts.getvdw12scale();
305 +    vdwScale_[2] = fopts.getvdw13scale();
306 +    vdwScale_[3] = fopts.getvdw14scale();
307 +    
308 +    electrostaticScale_[0] = 1.0;
309 +    electrostaticScale_[1] = fopts.getelectrostatic12scale();
310 +    electrostaticScale_[2] = fopts.getelectrostatic13scale();
311 +    electrostaticScale_[3] = fopts.getelectrostatic14scale();    
312 +    
313 +    fDecomp_->distributeInitialData();
314 +
315 +    initialized_ = true;
316 +
317 +  }
318 +
319 +  void ForceManager::calcForces() {
320 +    
321 +    if (!initialized_) initialize();
322 +
323 +    preCalculation();  
324 +    shortRangeInteractions();
325 +    longRangeInteractions();
326 +    postCalculation();    
327 +  }
328 +  
329    void ForceManager::preCalculation() {
330      SimInfo::MoleculeIterator mi;
331      Molecule* mol;
# Line 87 | Line 333 | namespace OpenMD {
333      Atom* atom;
334      Molecule::RigidBodyIterator rbIter;
335      RigidBody* rb;
336 +    Molecule::CutoffGroupIterator ci;
337 +    CutoffGroup* cg;
338      
339      // forces are zeroed here, before any are accumulated.
92 –    // NOTE: do not rezero the forces in Fortran.
340      
341      for (mol = info_->beginMolecule(mi); mol != NULL;
342           mol = info_->nextMolecule(mi)) {
# Line 102 | Line 349 | namespace OpenMD {
349             rb = mol->nextRigidBody(rbIter)) {
350          rb->zeroForcesAndTorques();
351        }        
352 <          
352 >
353 >      if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
354 >        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
355 >            cg = mol->nextCutoffGroup(ci)) {
356 >          //calculate the center of mass of cutoff group
357 >          cg->updateCOM();
358 >        }
359 >      }      
360      }
361 <    
361 >  
362      // Zero out the stress tensor
363      tau *= 0.0;
364      
365    }
366    
367 <  void ForceManager::calcShortRangeInteraction() {
367 >  void ForceManager::shortRangeInteractions() {
368      Molecule* mol;
369      RigidBody* rb;
370      Bond* bond;
# Line 152 | Line 406 | namespace OpenMD {
406          RealType currBendPot = bend->getPotential();          
407          
408          bendPotential += bend->getPotential();
409 <        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
409 >        map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
410          if (i == bendDataSets.end()) {
411            BendDataSet dataSet;
412            dataSet.prev.angle = dataSet.curr.angle = angle;
413            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
414            dataSet.deltaV = 0.0;
415 <          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
415 >          bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
416          }else {
417            i->second.prev.angle = i->second.curr.angle;
418            i->second.prev.potential = i->second.curr.potential;
# Line 175 | Line 429 | namespace OpenMD {
429          torsion->calcForce(angle);
430          RealType currTorsionPot = torsion->getPotential();
431          torsionPotential += torsion->getPotential();
432 <        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
432 >        map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
433          if (i == torsionDataSets.end()) {
434            TorsionDataSet dataSet;
435            dataSet.prev.angle = dataSet.curr.angle = angle;
436            dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
437            dataSet.deltaV = 0.0;
438 <          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
438 >          torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
439          }else {
440            i->second.prev.angle = i->second.curr.angle;
441            i->second.prev.potential = i->second.curr.potential;
# Line 191 | Line 445 | namespace OpenMD {
445                                     i->second.prev.potential);
446          }      
447        }      
448 <
448 >      
449        for (inversion = mol->beginInversion(inversionIter);
450             inversion != NULL;
451             inversion = mol->nextInversion(inversionIter)) {
# Line 199 | Line 453 | namespace OpenMD {
453          inversion->calcForce(angle);
454          RealType currInversionPot = inversion->getPotential();
455          inversionPotential += inversion->getPotential();
456 <        std::map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
456 >        map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
457          if (i == inversionDataSets.end()) {
458            InversionDataSet dataSet;
459            dataSet.prev.angle = dataSet.curr.angle = angle;
460            dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
461            dataSet.deltaV = 0.0;
462 <          inversionDataSets.insert(std::map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
462 >          inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
463          }else {
464            i->second.prev.angle = i->second.curr.angle;
465            i->second.prev.potential = i->second.curr.potential;
# Line 224 | Line 478 | namespace OpenMD {
478      curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
479      curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
480      curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
481 <    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
228 <    
481 >    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
482    }
483    
484 <  void ForceManager::calcLongRangeInteraction() {
232 <    Snapshot* curSnapshot;
233 <    DataStorage* config;
234 <    RealType* frc;
235 <    RealType* pos;
236 <    RealType* trq;
237 <    RealType* A;
238 <    RealType* electroFrame;
239 <    RealType* rc;
240 <    RealType* particlePot;
241 <    
242 <    //get current snapshot from SimInfo
243 <    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
244 <    
245 <    //get array pointers
246 <    config = &(curSnapshot->atomData);
247 <    frc = config->getArrayPointer(DataStorage::dslForce);
248 <    pos = config->getArrayPointer(DataStorage::dslPosition);
249 <    trq = config->getArrayPointer(DataStorage::dslTorque);
250 <    A   = config->getArrayPointer(DataStorage::dslAmat);
251 <    electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
252 <    particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
484 >  void ForceManager::longRangeInteractions() {
485  
486 +    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
487 +    DataStorage* config = &(curSnapshot->atomData);
488 +    DataStorage* cgConfig = &(curSnapshot->cgData);
489 +
490      //calculate the center of mass of cutoff group
491 +
492      SimInfo::MoleculeIterator mi;
493      Molecule* mol;
494      Molecule::CutoffGroupIterator ci;
495      CutoffGroup* cg;
496 <    Vector3d com;
497 <    std::vector<Vector3d> rcGroup;
261 <    
262 <    if(info_->getNCutoffGroups() > 0){
263 <      
496 >
497 >    if(info_->getNCutoffGroups() > 0){      
498        for (mol = info_->beginMolecule(mi); mol != NULL;
499             mol = info_->nextMolecule(mi)) {
500          for(cg = mol->beginCutoffGroup(ci); cg != NULL;
501              cg = mol->nextCutoffGroup(ci)) {
502 <          cg->getCOM(com);
269 <          rcGroup.push_back(com);
502 >          cg->updateCOM();
503          }
504 <      }// end for (mol)
272 <      
273 <      rc = rcGroup[0].getArrayPointer();
504 >      }      
505      } else {
506        // center of mass of the group is the same as position of the atom  
507        // if cutoff group does not exist
508 <      rc = pos;
508 >      cgConfig->position = config->position;
509      }
510 +
511 +    fDecomp_->zeroWorkArrays();
512 +    fDecomp_->distributeData();
513      
514 <    //initialize data before passing to fortran
515 <    RealType longRangePotential[LR_POT_TYPES];
516 <    RealType lrPot = 0.0;
517 <    int isError = 0;
514 >    int cg1, cg2, atom1, atom2, topoDist;
515 >    Vector3d d_grp, dag, d;
516 >    RealType rgrpsq, rgrp, r2, r;
517 >    RealType electroMult, vdwMult;
518 >    RealType vij;
519 >    Vector3d fij, fg, f1;
520 >    tuple3<RealType, RealType, RealType> cuts;
521 >    RealType rCutSq;
522 >    bool in_switching_region;
523 >    RealType sw, dswdr, swderiv;
524 >    vector<int> atomListColumn, atomListRow, atomListLocal;
525 >    InteractionData idat;
526 >    SelfData sdat;
527 >    RealType mf;
528 >    RealType lrPot;
529 >    RealType vpair;
530 >    potVec longRangePotential(0.0);
531 >    potVec workPot(0.0);
532  
533 <    for (int i=0; i<LR_POT_TYPES;i++){
534 <      longRangePotential[i]=0.0; //Initialize array
533 >    int loopStart, loopEnd;
534 >
535 >    idat.vdwMult = &vdwMult;
536 >    idat.electroMult = &electroMult;
537 >    idat.pot = &workPot;
538 >    sdat.pot = fDecomp_->getEmbeddingPotential();
539 >    idat.vpair = &vpair;
540 >    idat.f1 = &f1;
541 >    idat.sw = &sw;
542 >    idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
543 >    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
544 >    
545 >    loopEnd = PAIR_LOOP;
546 >    if (info_->requiresPrepair() ) {
547 >      loopStart = PREPAIR_LOOP;
548 >    } else {
549 >      loopStart = PAIR_LOOP;
550      }
551 +  
552 +    for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
553      
554 <    doForceLoop(pos,
555 <                rc,
556 <                A,
557 <                electroFrame,
558 <                frc,
559 <                trq,
560 <                tau.getArrayPointer(),
561 <                longRangePotential,
562 <                particlePot,
563 <                &isError );
554 >      if (iLoop == loopStart) {
555 >        bool update_nlist = fDecomp_->checkNeighborList();
556 >        if (update_nlist)
557 >          neighborList = fDecomp_->buildNeighborList();
558 >      }      
559 >        
560 >      for (vector<pair<int, int> >::iterator it = neighborList.begin();
561 >             it != neighborList.end(); ++it) {
562 >                
563 >        cg1 = (*it).first;
564 >        cg2 = (*it).second;
565 >        
566 >        cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
567 >
568 >        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
569 >        curSnapshot->wrapVector(d_grp);        
570 >        rgrpsq = d_grp.lengthSquare();
571 >
572 >        rCutSq = cuts.second;
573 >
574 >        if (rgrpsq < rCutSq) {
575 >          idat.rcut = &cuts.first;
576 >          if (iLoop == PAIR_LOOP) {
577 >            vij *= 0.0;
578 >            fij = V3Zero;
579 >          }
580 >          
581 >          in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
582 >                                                     rgrp);
583 >              
584 >          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
585 >          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
586 >
587 >          for (vector<int>::iterator ia = atomListRow.begin();
588 >               ia != atomListRow.end(); ++ia) {            
589 >            atom1 = (*ia);
590 >            
591 >            for (vector<int>::iterator jb = atomListColumn.begin();
592 >                 jb != atomListColumn.end(); ++jb) {              
593 >              atom2 = (*jb);
594 >            
595 >              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
596 >                vpair = 0.0;
597 >                workPot = 0.0;
598 >                f1 = V3Zero;
599 >
600 >                fDecomp_->fillInteractionData(idat, atom1, atom2);
601 >                
602 >                topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
603 >                vdwMult = vdwScale_[topoDist];
604 >                electroMult = electrostaticScale_[topoDist];
605 >
606 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
607 >                  idat.d = &d_grp;
608 >                  idat.r2 = &rgrpsq;
609 >                } else {
610 >                  d = fDecomp_->getInteratomicVector(atom1, atom2);
611 >                  curSnapshot->wrapVector( d );
612 >                  r2 = d.lengthSquare();
613 >                  idat.d = &d;
614 >                  idat.r2 = &r2;
615 >                }
616 >                
617 >                r = sqrt( *(idat.r2) );
618 >                idat.rij = &r;
619 >              
620 >                if (iLoop == PREPAIR_LOOP) {
621 >                  interactionMan_->doPrePair(idat);
622 >                } else {
623 >                  interactionMan_->doPair(idat);
624 >                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
625 >                  vij += vpair;
626 >                  fij += f1;
627 >                  tau -= outProduct( *(idat.d), f1);
628 >                }
629 >              }
630 >            }
631 >          }
632 >
633 >          if (iLoop == PAIR_LOOP) {
634 >            if (in_switching_region) {
635 >              swderiv = vij * dswdr / rgrp;
636 >              fg = swderiv * d_grp;
637 >
638 >              fij += fg;
639 >
640 >              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
641 >                tau -= outProduct( *(idat.d), fg);
642 >              }
643 >          
644 >              for (vector<int>::iterator ia = atomListRow.begin();
645 >                   ia != atomListRow.end(); ++ia) {            
646 >                atom1 = (*ia);                
647 >                mf = fDecomp_->getMassFactorRow(atom1);
648 >                // fg is the force on atom ia due to cutoff group's
649 >                // presence in switching region
650 >                fg = swderiv * d_grp * mf;
651 >                fDecomp_->addForceToAtomRow(atom1, fg);
652 >
653 >                if (atomListRow.size() > 1) {
654 >                  if (info_->usesAtomicVirial()) {
655 >                    // find the distance between the atom
656 >                    // and the center of the cutoff group:
657 >                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
658 >                    tau -= outProduct(dag, fg);
659 >                  }
660 >                }
661 >              }
662 >              for (vector<int>::iterator jb = atomListColumn.begin();
663 >                   jb != atomListColumn.end(); ++jb) {              
664 >                atom2 = (*jb);
665 >                mf = fDecomp_->getMassFactorColumn(atom2);
666 >                // fg is the force on atom jb due to cutoff group's
667 >                // presence in switching region
668 >                fg = -swderiv * d_grp * mf;
669 >                fDecomp_->addForceToAtomColumn(atom2, fg);
670 >
671 >                if (atomListColumn.size() > 1) {
672 >                  if (info_->usesAtomicVirial()) {
673 >                    // find the distance between the atom
674 >                    // and the center of the cutoff group:
675 >                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
676 >                    tau -= outProduct(dag, fg);
677 >                  }
678 >                }
679 >              }
680 >            }
681 >            //if (!SIM_uses_AtomicVirial) {
682 >            //  tau -= outProduct(d_grp, fij);
683 >            //}
684 >          }
685 >        }
686 >      }
687 >
688 >      if (iLoop == PREPAIR_LOOP) {
689 >        if (info_->requiresPrepair()) {            
690 >          fDecomp_->collectIntermediateData();
691 >
692 >          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
693 >            fDecomp_->fillSelfData(sdat, atom1);
694 >            interactionMan_->doPreForce(sdat);
695 >          }
696 >          
697 >          
698 >          fDecomp_->distributeIntermediateData();        
699 >        }
700 >      }
701 >
702 >    }
703      
704 <    if( isError ){
705 <      sprintf( painCave.errMsg,
706 <               "Error returned from the fortran force calculation.\n" );
707 <      painCave.isFatal = 1;
708 <      simError();
704 >    fDecomp_->collectData();
705 >    
706 >    if ( info_->requiresSkipCorrection() ) {
707 >      
708 >      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
709 >
710 >        vector<int> skipList = fDecomp_->getSkipsForAtom( atom1 );
711 >        
712 >        for (vector<int>::iterator jb = skipList.begin();
713 >             jb != skipList.end(); ++jb) {        
714 >    
715 >          atom2 = (*jb);
716 >          fDecomp_->fillSkipData(idat, atom1, atom2);
717 >          interactionMan_->doSkipCorrection(idat);
718 >          fDecomp_->unpackSkipData(idat, atom1, atom2);
719 >
720 >        }
721 >      }
722      }
723 <    for (int i=0; i<LR_POT_TYPES;i++){
724 <      lrPot += longRangePotential[i]; //Quick hack
723 >    
724 >    if (info_->requiresSelfCorrection()) {
725 >
726 >      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
727 >        fDecomp_->fillSelfData(sdat, atom1);
728 >        interactionMan_->doSelfCorrection(sdat);
729 >      }
730 >
731      }
732 <        
732 >
733 >    longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
734 >      *(fDecomp_->getPairwisePotential());
735 >
736 >    lrPot = longRangePotential.sum();
737 >
738      //store the tau and long range potential    
739      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
740 <    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
741 <    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
740 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
741 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
742    }
743  
744    

Diff Legend

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