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

Comparing:
trunk/src/brains/ForceManager.cpp (file contents), Revision 1126 by gezelter, Fri Apr 6 21:53:43 2007 UTC vs.
branches/development/src/brains/ForceManager.cpp (file contents), Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   /**
# Line 49 | Line 49
49  
50   #include "brains/ForceManager.hpp"
51   #include "primitives/Molecule.hpp"
52 < #include "UseTheForce/doForces_interface.h"
53 < #define __C
54 < #include "UseTheForce/DarkSide/fInteractionMap.h"
52 > #define __OPENMD_C
53   #include "utils/simError.h"
54 + #include "primitives/Bond.hpp"
55   #include "primitives/Bend.hpp"
56 < #include "primitives/Bend.hpp"
57 < namespace oopse {
56 > #include "primitives/Torsion.hpp"
57 > #include "primitives/Inversion.hpp"
58 > #include "nonbonded/NonBondedInteraction.hpp"
59 > #include "parallel/ForceMatrixDecomposition.hpp"
60  
61 <  void ForceManager::calcForces(bool needPotential, bool needStress) {
61 > using namespace std;
62 > namespace OpenMD {
63 >  
64 >  ForceManager::ForceManager(SimInfo * info) : info_(info) {
65 >
66 >    fDecomp_ = new ForceMatrixDecomposition(info_);
67 >  }
68 >  
69 >  void ForceManager::calcForces() {
70      
71 <    if (!info_->isFortranInitialized()) {
71 >    if (!info_->isTopologyDone()) {
72        info_->update();
73 +      interactionMan_->setSimInfo(info_);
74 +      interactionMan_->initialize();
75 +      swfun_ = interactionMan_->getSwitchingFunction();
76 +      info_->prepareTopology();
77 +      fDecomp_->distributeInitialData();
78      }
79      
80 <    preCalculation();
80 >    preCalculation();  
81 >    shortRangeInteractions();
82 >    longRangeInteractions();
83 >    postCalculation();
84      
68    calcShortRangeInteraction();
69
70    calcLongRangeInteraction(needPotential, needStress);
71
72    postCalculation(needStress);
73    
85    }
86    
87    void ForceManager::preCalculation() {
# Line 80 | Line 91 | namespace oopse {
91      Atom* atom;
92      Molecule::RigidBodyIterator rbIter;
93      RigidBody* rb;
94 +    Molecule::CutoffGroupIterator ci;
95 +    CutoffGroup* cg;
96      
97      // forces are zeroed here, before any are accumulated.
98 <    // NOTE: do not rezero the forces in Fortran.
86 <
98 >    
99      for (mol = info_->beginMolecule(mi); mol != NULL;
100           mol = info_->nextMolecule(mi)) {
101        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
102          atom->zeroForcesAndTorques();
103        }
104 <      
104 >          
105        //change the positions of atoms which belong to the rigidbodies
106        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
107             rb = mol->nextRigidBody(rbIter)) {
108          rb->zeroForcesAndTorques();
109        }        
110 +
111 +      if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
112 +        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
113 +            cg = mol->nextCutoffGroup(ci)) {
114 +          //calculate the center of mass of cutoff group
115 +          cg->updateCOM();
116 +        }
117 +      }      
118      }
119 <    
119 >  
120      // Zero out the stress tensor
121      tau *= 0.0;
122      
123    }
124    
125 <  void ForceManager::calcShortRangeInteraction() {
125 >  void ForceManager::shortRangeInteractions() {
126      Molecule* mol;
127      RigidBody* rb;
128      Bond* bond;
129      Bend* bend;
130      Torsion* torsion;
131 +    Inversion* inversion;
132      SimInfo::MoleculeIterator mi;
133      Molecule::RigidBodyIterator rbIter;
134      Molecule::BondIterator bondIter;;
135      Molecule::BendIterator  bendIter;
136      Molecule::TorsionIterator  torsionIter;
137 +    Molecule::InversionIterator  inversionIter;
138      RealType bondPotential = 0.0;
139      RealType bendPotential = 0.0;
140      RealType torsionPotential = 0.0;
141 +    RealType inversionPotential = 0.0;
142  
143      //calculate short range interactions    
144      for (mol = info_->beginMolecule(mi); mol != NULL;
# Line 139 | Line 162 | namespace oopse {
162          RealType angle;
163          bend->calcForce(angle);
164          RealType currBendPot = bend->getPotential();          
165 +        
166          bendPotential += bend->getPotential();
167 <        std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
167 >        map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
168          if (i == bendDataSets.end()) {
169            BendDataSet dataSet;
170            dataSet.prev.angle = dataSet.curr.angle = angle;
171            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
172            dataSet.deltaV = 0.0;
173 <          bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
173 >          bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
174          }else {
175            i->second.prev.angle = i->second.curr.angle;
176            i->second.prev.potential = i->second.curr.potential;
# Line 163 | Line 187 | namespace oopse {
187          torsion->calcForce(angle);
188          RealType currTorsionPot = torsion->getPotential();
189          torsionPotential += torsion->getPotential();
190 <        std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
190 >        map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
191          if (i == torsionDataSets.end()) {
192            TorsionDataSet dataSet;
193            dataSet.prev.angle = dataSet.curr.angle = angle;
194            dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
195            dataSet.deltaV = 0.0;
196 <          torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
196 >          torsionDataSets.insert(map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
197          }else {
198            i->second.prev.angle = i->second.curr.angle;
199            i->second.prev.potential = i->second.curr.potential;
# Line 179 | Line 203 | namespace oopse {
203                                     i->second.prev.potential);
204          }      
205        }      
206 +      
207 +      for (inversion = mol->beginInversion(inversionIter);
208 +           inversion != NULL;
209 +           inversion = mol->nextInversion(inversionIter)) {
210 +        RealType angle;
211 +        inversion->calcForce(angle);
212 +        RealType currInversionPot = inversion->getPotential();
213 +        inversionPotential += inversion->getPotential();
214 +        map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
215 +        if (i == inversionDataSets.end()) {
216 +          InversionDataSet dataSet;
217 +          dataSet.prev.angle = dataSet.curr.angle = angle;
218 +          dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
219 +          dataSet.deltaV = 0.0;
220 +          inversionDataSets.insert(map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
221 +        }else {
222 +          i->second.prev.angle = i->second.curr.angle;
223 +          i->second.prev.potential = i->second.curr.potential;
224 +          i->second.curr.angle = angle;
225 +          i->second.curr.potential = currInversionPot;
226 +          i->second.deltaV =  fabs(i->second.curr.potential -  
227 +                                   i->second.prev.potential);
228 +        }      
229 +      }      
230      }
231      
232      RealType  shortRangePotential = bondPotential + bendPotential +
233 <      torsionPotential;    
233 >      torsionPotential +  inversionPotential;    
234      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
235      curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
236      curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
237      curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
238      curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
239 <    
239 >    curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;    
240    }
241    
242 <  void ForceManager::calcLongRangeInteraction(bool needPotential,
195 <                                              bool needStress) {
196 <    Snapshot* curSnapshot;
197 <    DataStorage* config;
198 <    RealType* frc;
199 <    RealType* pos;
200 <    RealType* trq;
201 <    RealType* A;
202 <    RealType* electroFrame;
203 <    RealType* rc;
204 <    
205 <    //get current snapshot from SimInfo
206 <    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
207 <    
208 <    //get array pointers
209 <    config = &(curSnapshot->atomData);
210 <    frc = config->getArrayPointer(DataStorage::dslForce);
211 <    pos = config->getArrayPointer(DataStorage::dslPosition);
212 <    trq = config->getArrayPointer(DataStorage::dslTorque);
213 <    A   = config->getArrayPointer(DataStorage::dslAmat);
214 <    electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
242 >  void ForceManager::longRangeInteractions() {
243  
244 <    //calculate the center of mass of cutoff group
245 <    SimInfo::MoleculeIterator mi;
246 <    Molecule* mol;
247 <    Molecule::CutoffGroupIterator ci;
248 <    CutoffGroup* cg;
249 <    Vector3d com;
250 <    std::vector<Vector3d> rcGroup;
251 <    
252 <    if(info_->getNCutoffGroups() > 0){
253 <      
254 <      for (mol = info_->beginMolecule(mi); mol != NULL;
255 <           mol = info_->nextMolecule(mi)) {
256 <        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
257 <            cg = mol->nextCutoffGroup(ci)) {
230 <          cg->getCOM(com);
231 <          rcGroup.push_back(com);
232 <        }
233 <      }// end for (mol)
234 <      
235 <      rc = rcGroup[0].getArrayPointer();
244 >    // some of this initial stuff will go away:
245 >    Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
246 >    DataStorage* config = &(curSnapshot->atomData);
247 >    DataStorage* cgConfig = &(curSnapshot->cgData);
248 >    RealType* frc = config->getArrayPointer(DataStorage::dslForce);
249 >    RealType* pos = config->getArrayPointer(DataStorage::dslPosition);
250 >    RealType* trq = config->getArrayPointer(DataStorage::dslTorque);
251 >    RealType* A = config->getArrayPointer(DataStorage::dslAmat);
252 >    RealType* electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
253 >    RealType* particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
254 >    RealType* rc;    
255 >
256 >    if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
257 >      rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
258      } else {
259        // center of mass of the group is the same as position of the atom  
260        // if cutoff group does not exist
# Line 240 | Line 262 | namespace oopse {
262      }
263      
264      //initialize data before passing to fortran
265 <    RealType longRangePotential[LR_POT_TYPES];
265 >    RealType longRangePotential[N_INTERACTION_FAMILIES];
266      RealType lrPot = 0.0;
245    Vector3d totalDipole;
246    short int passedCalcPot = needPotential;
247    short int passedCalcStress = needStress;
267      int isError = 0;
268  
269 <    for (int i=0; i<LR_POT_TYPES;i++){
269 >    // dangerous to iterate over enums, but we'll live on the edge:
270 >    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
271        longRangePotential[i]=0.0; //Initialize array
272      }
253    
254    doForceLoop( pos,
255                 rc,
256                 A,
257                 electroFrame,
258                 frc,
259                 trq,
260                 tau.getArrayPointer(),
261                 longRangePotential,
262                 &passedCalcPot,
263                 &passedCalcStress,
264                 &isError );
273  
274 <    if( isError ){
275 <      sprintf( painCave.errMsg,
276 <               "Error returned from the fortran force calculation.\n" );
277 <      painCave.isFatal = 1;
278 <      simError();
274 >    // new stuff starts here:
275 >
276 >    fDecomp_->distributeData();
277 >
278 >    int cg1, cg2, atom1, atom2;
279 >    Vector3d d_grp, dag;
280 >    RealType rgrpsq, rgrp;
281 >    RealType vij;
282 >    Vector3d fij, fg;
283 >    pair<int, int> gtypes;
284 >    RealType rCutSq;
285 >    bool in_switching_region;
286 >    RealType sw, dswdr, swderiv;
287 >    vector<int> atomListColumn, atomListRow, atomListLocal;
288 >    InteractionData idat;
289 >    SelfData sdat;
290 >    RealType mf;
291 >
292 >    int loopStart, loopEnd;
293 >
294 >    loopEnd = PAIR_LOOP;
295 >    if (info_->requiresPrepair() ) {
296 >      loopStart = PREPAIR_LOOP;
297 >    } else {
298 >      loopStart = PAIR_LOOP;
299      }
300 <    for (int i=0; i<LR_POT_TYPES;i++){
301 <      lrPot += longRangePotential[i]; //Quick hack
300 >
301 >    for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
302 >      
303 >      if (iLoop == loopStart) {
304 >        bool update_nlist = fDecomp_->checkNeighborList();
305 >        if (update_nlist)
306 >          neighborList = fDecomp_->buildNeighborList();
307 >      }
308 >
309 >      for (vector<pair<int, int> >::iterator it = neighborList.begin();
310 >             it != neighborList.end(); ++it) {
311 >        
312 >        cg1 = (*it).first;
313 >        cg2 = (*it).second;
314 >
315 >        gtypes = fDecomp_->getGroupTypes(cg1, cg2);
316 >        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
317 >        curSnapshot->wrapVector(d_grp);        
318 >        rgrpsq = d_grp.lengthSquare();
319 >        rCutSq = groupCutoffMap[gtypes].first;
320 >
321 >        if (rgrpsq < rCutSq) {
322 >          *(idat.rcut) = groupCutoffMap[gtypes].second;
323 >          if (iLoop == PAIR_LOOP) {
324 >            vij *= 0.0;
325 >            fij = V3Zero;
326 >          }
327 >          
328 >          in_switching_region = swfun_->getSwitch(rgrpsq, *(idat.sw), dswdr,
329 >                                                  rgrp);              
330 >          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
331 >          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
332 >
333 >          for (vector<int>::iterator ia = atomListRow.begin();
334 >               ia != atomListRow.end(); ++ia) {            
335 >            atom1 = (*ia);
336 >            
337 >            for (vector<int>::iterator jb = atomListColumn.begin();
338 >                 jb != atomListColumn.end(); ++jb) {              
339 >              atom2 = (*jb);
340 >              
341 >              if (!fDecomp_->skipAtomPair(atom1, atom2)) {
342 >                
343 >                idat = fDecomp_->fillInteractionData(atom1, atom2);
344 >
345 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
346 >                  *(idat.d) = d_grp;
347 >                  *(idat.r2) = rgrpsq;
348 >                } else {
349 >                  *(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2);
350 >                  curSnapshot->wrapVector( *(idat.d) );
351 >                  *(idat.r2) = idat.d->lengthSquare();
352 >                }
353 >                
354 >                *(idat.rij) = sqrt( *(idat.r2) );
355 >              
356 >                if (iLoop == PREPAIR_LOOP) {
357 >                  interactionMan_->doPrePair(idat);
358 >                } else {
359 >                  interactionMan_->doPair(idat);
360 >                  vij += *(idat.vpair);
361 >                  fij += *(idat.f1);
362 >                  tau -= outProduct( *(idat.d), *(idat.f1));
363 >                }
364 >              }
365 >            }
366 >          }
367 >
368 >          if (iLoop == PAIR_LOOP) {
369 >            if (in_switching_region) {
370 >              swderiv = vij * dswdr / rgrp;
371 >              fg = swderiv * d_grp;
372 >
373 >              fij += fg;
374 >
375 >              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
376 >                tau -= outProduct( *(idat.d), fg);
377 >              }
378 >          
379 >              for (vector<int>::iterator ia = atomListRow.begin();
380 >                   ia != atomListRow.end(); ++ia) {            
381 >                atom1 = (*ia);                
382 >                mf = fDecomp_->getMassFactorRow(atom1);
383 >                // fg is the force on atom ia due to cutoff group's
384 >                // presence in switching region
385 >                fg = swderiv * d_grp * mf;
386 >                fDecomp_->addForceToAtomRow(atom1, fg);
387 >
388 >                if (atomListRow.size() > 1) {
389 >                  if (info_->usesAtomicVirial()) {
390 >                    // find the distance between the atom
391 >                    // and the center of the cutoff group:
392 >                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
393 >                    tau -= outProduct(dag, fg);
394 >                  }
395 >                }
396 >              }
397 >              for (vector<int>::iterator jb = atomListColumn.begin();
398 >                   jb != atomListColumn.end(); ++jb) {              
399 >                atom2 = (*jb);
400 >                mf = fDecomp_->getMassFactorColumn(atom2);
401 >                // fg is the force on atom jb due to cutoff group's
402 >                // presence in switching region
403 >                fg = -swderiv * d_grp * mf;
404 >                fDecomp_->addForceToAtomColumn(atom2, fg);
405 >
406 >                if (atomListColumn.size() > 1) {
407 >                  if (info_->usesAtomicVirial()) {
408 >                    // find the distance between the atom
409 >                    // and the center of the cutoff group:
410 >                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
411 >                    tau -= outProduct(dag, fg);
412 >                  }
413 >                }
414 >              }
415 >            }
416 >            //if (!SIM_uses_AtomicVirial) {
417 >            //  tau -= outProduct(d_grp, fij);
418 >            //}
419 >          }
420 >        }
421 >      }
422 >
423 >      if (iLoop == PREPAIR_LOOP) {
424 >        if (info_->requiresPrepair()) {            
425 >          fDecomp_->collectIntermediateData();
426 >
427 >          for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
428 >            sdat = fDecomp_->fillSelfData(atom1);
429 >            interactionMan_->doPreForce(sdat);
430 >          }
431 >
432 >          fDecomp_->distributeIntermediateData();        
433 >        }
434 >      }
435 >
436      }
437      
438 <    // grab the simulation box dipole moment if specified
439 <    if (info_->getCalcBoxDipole()){
440 <      getAccumulatedBoxDipole(totalDipole.getArrayPointer());
438 >    fDecomp_->collectData();
439 >    
440 >    if ( info_->requiresSkipCorrection() ) {
441        
442 <      curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0);
443 <      curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1);
444 <      curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2);
442 >      for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
443 >
444 >        vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
445 >        
446 >        for (vector<int>::iterator jb = skipList.begin();
447 >             jb != skipList.end(); ++jb) {        
448 >    
449 >          atom2 = (*jb);
450 >          idat = fDecomp_->fillSkipData(atom1, atom2);
451 >          interactionMan_->doSkipCorrection(idat);
452 >
453 >        }
454 >      }
455      }
456      
457 +    if (info_->requiresSelfCorrection()) {
458 +
459 +      for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {          
460 +        sdat = fDecomp_->fillSelfData(atom1);
461 +        interactionMan_->doSelfCorrection(sdat);
462 +      }
463 +
464 +    }
465 +
466 +    // dangerous to iterate over enums, but we'll live on the edge:
467 +    for (int i = NO_FAMILY; i != N_INTERACTION_FAMILIES; ++i){
468 +      lrPot += longRangePotential[i]; //Quick hack
469 +    }
470 +        
471      //store the tau and long range potential    
472      curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
473 <    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
474 <    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
473 >    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
474 >    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
475    }
476  
477    
478 <  void ForceManager::postCalculation(bool needStress) {
478 >  void ForceManager::postCalculation() {
479      SimInfo::MoleculeIterator mi;
480      Molecule* mol;
481      Molecule::RigidBodyIterator rbIter;
# Line 302 | Line 488 | namespace oopse {
488           mol = info_->nextMolecule(mi)) {
489        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
490             rb = mol->nextRigidBody(rbIter)) {
491 <        if (needStress) {          
492 <          Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
307 <          tau += rbTau;
308 <        } else{
309 <          rb->calcForcesAndTorques();
310 <        }
491 >        Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
492 >        tau += rbTau;
493        }
494      }
495 <
314 <    if (needStress) {
495 >    
496   #ifdef IS_MPI
497 <      Mat3x3d tmpTau(tau);
498 <      MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
499 <                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
497 >    Mat3x3d tmpTau(tau);
498 >    MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
499 >                  9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
500   #endif
501 <      curSnapshot->statData.setTau(tau);
321 <    }
501 >    curSnapshot->statData.setTau(tau);
502    }
503  
504 < } //end namespace oopse
504 > } //end namespace OpenMD

Comparing:
trunk/src/brains/ForceManager.cpp (property svn:keywords), Revision 1126 by gezelter, Fri Apr 6 21:53:43 2007 UTC vs.
branches/development/src/brains/ForceManager.cpp (property svn:keywords), Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines