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

Comparing trunk/src/brains/ForceManager.cpp (file contents):
Revision 1925 by gezelter, Wed Aug 7 15:24:16 2013 UTC vs.
Revision 2060 by gezelter, Tue Mar 3 16:19:47 2015 UTC

# Line 57 | Line 57
57   #include "primitives/Torsion.hpp"
58   #include "primitives/Inversion.hpp"
59   #include "nonbonded/NonBondedInteraction.hpp"
60 < #include "perturbations/ElectricField.hpp"
60 > #include "perturbations/UniformField.hpp"
61 > #include "perturbations/UniformGradient.hpp"
62   #include "parallel/ForceMatrixDecomposition.hpp"
63  
64   #include <cstdio>
# Line 87 | Line 88 | namespace OpenMD {
88    /**
89     * setupCutoffs
90     *
91 <   * Sets the values of cutoffRadius, switchingRadius, cutoffMethod,
91 <   * and cutoffPolicy
91 >   * Sets the values of cutoffRadius, switchingRadius, and cutoffMethod
92     *
93     * cutoffRadius : realType
94     *  If the cutoffRadius was explicitly set, use that value.
# Line 103 | Line 103 | namespace OpenMD {
103     *      If cutoffMethod was explicitly set, use that choice.
104     *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
105     *
106 –   * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
107 –   *      If cutoffPolicy was explicitly set, use that choice.
108 –   *      If cutoffPolicy was not explicitly set, use TRADITIONAL
109 –   *
106     * switchingRadius : realType
107     *  If the cutoffMethod was set to SWITCHED:
108     *      If the switchingRadius was explicitly set, use that value
# Line 163 | Line 159 | namespace OpenMD {
159        }
160      }
161  
162 <    fDecomp_->setUserCutoff(rCut_);
162 >    fDecomp_->setCutoffRadius(rCut_);
163      interactionMan_->setCutoffRadius(rCut_);
164 +    rCutSq_ = rCut_ * rCut_;
165  
166      map<string, CutoffMethod> stringToCutoffMethod;
167      stringToCutoffMethod["HARD"] = HARD;
# Line 273 | Line 270 | namespace OpenMD {
270            }
271          }
272        }
276 –    }
277 –
278 –    map<string, CutoffPolicy> stringToCutoffPolicy;
279 –    stringToCutoffPolicy["MIX"] = MIX;
280 –    stringToCutoffPolicy["MAX"] = MAX;
281 –    stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;    
282 –
283 –    string cutPolicy;
284 –    if (forceFieldOptions_.haveCutoffPolicy()){
285 –      cutPolicy = forceFieldOptions_.getCutoffPolicy();
286 –    }else if (simParams_->haveCutoffPolicy()) {
287 –      cutPolicy = simParams_->getCutoffPolicy();
273      }
289 –
290 –    if (!cutPolicy.empty()){
291 –      toUpper(cutPolicy);
292 –      map<string, CutoffPolicy>::iterator i;
293 –      i = stringToCutoffPolicy.find(cutPolicy);
294 –
295 –      if (i == stringToCutoffPolicy.end()) {
296 –        sprintf(painCave.errMsg,
297 –                "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
298 –                "\tShould be one of: "
299 –                "MIX, MAX, or TRADITIONAL\n",
300 –                cutPolicy.c_str());
301 –        painCave.isFatal = 1;
302 –        painCave.severity = OPENMD_ERROR;
303 –        simError();
304 –      } else {
305 –        cutoffPolicy_ = i->second;
306 –      }
307 –    } else {
308 –      sprintf(painCave.errMsg,
309 –              "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
310 –              "\tOpenMD will use TRADITIONAL.\n");
311 –      painCave.isFatal = 0;
312 –      painCave.severity = OPENMD_INFO;
313 –      simError();
314 –      cutoffPolicy_ = TRADITIONAL;        
315 –    }
316 –
317 –    fDecomp_->setCutoffPolicy(cutoffPolicy_);
274          
275      // create the switching function object:
276  
# Line 394 | Line 350 | namespace OpenMD {
350      switcher_->setSwitch(rSwitch_, rCut_);
351    }
352  
397 –
398 –
399 –  
353    void ForceManager::initialize() {
354  
355      if (!info_->isTopologyDone()) {
# Line 405 | Line 358 | namespace OpenMD {
358        interactionMan_->setSimInfo(info_);
359        interactionMan_->initialize();
360  
361 <      // We want to delay the cutoffs until after the interaction
362 <      // manager has set up the atom-atom interactions so that we can
363 <      // query them for suggested cutoff values
361 >      //! We want to delay the cutoffs until after the interaction
362 >      //! manager has set up the atom-atom interactions so that we can
363 >      //! query them for suggested cutoff values
364        setupCutoffs();
365  
366        info_->prepareTopology();      
# Line 417 | Line 370 | namespace OpenMD {
370        if (doHeatFlux_) doParticlePot_ = true;
371  
372        doElectricField_ = info_->getSimParams()->getOutputElectricField();
373 +      doSitePotential_ = info_->getSimParams()->getOutputSitePotential();
374    
375      }
376  
377      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
378      
379 <    // Force fields can set options on how to scale van der Waals and
380 <    // electrostatic interactions for atoms connected via bonds, bends
381 <    // and torsions in this case the topological distance between
382 <    // atoms is:
383 <    // 0 = topologically unconnected
384 <    // 1 = bonded together
385 <    // 2 = connected via a bend
386 <    // 3 = connected via a torsion
379 >    //! Force fields can set options on how to scale van der Waals and
380 >    //! electrostatic interactions for atoms connected via bonds, bends
381 >    //! and torsions in this case the topological distance between
382 >    //! atoms is:
383 >    //! 0 = topologically unconnected
384 >    //! 1 = bonded together
385 >    //! 2 = connected via a bend
386 >    //! 3 = connected via a torsion
387      
388      vdwScale_.reserve(4);
389      fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
# Line 447 | Line 401 | namespace OpenMD {
401      electrostaticScale_[2] = fopts.getelectrostatic13scale();
402      electrostaticScale_[3] = fopts.getelectrostatic14scale();    
403      
404 <    if (info_->getSimParams()->haveElectricField()) {
405 <      ElectricField* eField = new ElectricField(info_);
404 >    if (info_->getSimParams()->haveUniformField()) {
405 >      UniformField* eField = new UniformField(info_);
406        perturbations_.push_back(eField);
407      }
408 <
408 >    if (info_->getSimParams()->haveUniformGradientStrength() ||
409 >        info_->getSimParams()->haveUniformGradientDirection1() ||
410 >        info_->getSimParams()->haveUniformGradientDirection2() ) {
411 >      UniformGradient* eGrad = new UniformGradient(info_);
412 >      perturbations_.push_back(eGrad);
413 >    }
414 >    
415      usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions();
416      
417      fDecomp_->distributeInitialData();
# Line 637 | Line 597 | namespace OpenMD {
597      // Collect from all nodes.  This should eventually be moved into a
598      // SystemDecomposition, but this is a better place than in
599      // Thermo to do the collection.
600 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE,
601 <                              MPI::SUM);
602 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE,
603 <                              MPI::SUM);
604 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1,
605 <                              MPI::REALTYPE, MPI::SUM);
606 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1,
607 <                              MPI::REALTYPE, MPI::SUM);
600 >
601 >    MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE,
602 >                  MPI_SUM, MPI_COMM_WORLD);
603 >    MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE,
604 >                  MPI_SUM, MPI_COMM_WORLD);
605 >    MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1,
606 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
607 >    MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1,
608 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
609   #endif
610  
611      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 666 | Line 627 | namespace OpenMD {
627      DataStorage* config = &(curSnapshot->atomData);
628      DataStorage* cgConfig = &(curSnapshot->cgData);
629  
669 –
630      //calculate the center of mass of cutoff group
631  
632      SimInfo::MoleculeIterator mi;
633      Molecule* mol;
634      Molecule::CutoffGroupIterator ci;
635      CutoffGroup* cg;
636 <
637 <    if(info_->getNCutoffGroups() > 0){      
636 >    
637 >    if(info_->getNCutoffGroups() != info_->getNAtoms()){
638        for (mol = info_->beginMolecule(mi); mol != NULL;
639             mol = info_->nextMolecule(mi)) {
640          for(cg = mol->beginCutoffGroup(ci); cg != NULL;
# Line 698 | Line 658 | namespace OpenMD {
658      RealType electroMult, vdwMult;
659      RealType vij;
660      Vector3d fij, fg, f1;
701 –    tuple3<RealType, RealType, RealType> cuts;
702 –    RealType rCut, rCutSq, rListSq;
661      bool in_switching_region;
662      RealType sw, dswdr, swderiv;
663      vector<int> atomListColumn, atomListRow;
# Line 715 | Line 673 | namespace OpenMD {
673      potVec exPot(0.0);
674      Vector3d eField1(0.0);
675      Vector3d eField2(0.0);
676 +    RealType sPot1(0.0);
677 +    RealType sPot2(0.0);
678 +    bool newAtom1;
679 +                  
680      vector<int>::iterator ia, jb;
681  
682      int loopStart, loopEnd;
683      
684 <    idat.rcut = &rCut;
684 >    idat.rcut = &rCut_;
685      idat.vdwMult = &vdwMult;
686      idat.electroMult = &electroMult;
687      idat.pot = &workPot;
# Line 730 | Line 692 | namespace OpenMD {
692      idat.dVdFQ1 = &dVdFQ1;
693      idat.dVdFQ2 = &dVdFQ2;
694      idat.eField1 = &eField1;
695 <    idat.eField2 = &eField2;  
695 >    idat.eField2 = &eField2;
696 >    idat.sPot1 = &sPot1;
697 >    idat.sPot2 = &sPot2;
698      idat.f1 = &f1;
699      idat.sw = &sw;
700      idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
701 <    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE || cutoffMethod_ == TAYLOR_SHIFTED) ? true : false;
701 >    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE ||
702 >                         cutoffMethod_ == TAYLOR_SHIFTED) ? true : false;
703      idat.doParticlePot = doParticlePot_;
704      idat.doElectricField = doElectricField_;
705 +    idat.doSitePotential = doSitePotential_;
706      sdat.doParticlePot = doParticlePot_;
707      
708      loopEnd = PAIR_LOOP;
# Line 752 | Line 718 | namespace OpenMD {
718          if (update_nlist) {
719            if (!usePeriodicBoundaryConditions_)
720              Mat3x3d bbox = thermo->getBoundingBox();
721 <          fDecomp_->buildNeighborList(neighborList_);
721 >          fDecomp_->buildNeighborList(neighborList_, point_);
722          }
723        }
724  
725 <      for (vector<pair<int, int> >::iterator it = neighborList_.begin();
760 <             it != neighborList_.end(); ++it) {
761 <                
762 <        cg1 = (*it).first;
763 <        cg2 = (*it).second;
725 >      for (cg1 = 0; cg1 < point_.size() - 1; cg1++) {
726          
727 <        fDecomp_->getGroupCutoffs(cg1, cg2, rCut, rCutSq, rListSq);
727 >        atomListRow = fDecomp_->getAtomsInGroupRow(cg1);        
728 >        newAtom1 = true;
729 >        
730 >        for (int m2 = point_[cg1]; m2 < point_[cg1+1]; m2++) {
731  
732 <        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
768 <
769 <        // already wrapped in the getIntergroupVector call:
770 <        // curSnapshot->wrapVector(d_grp);        
771 <        rgrpsq = d_grp.lengthSquare();
772 <
773 <        if (rgrpsq < rCutSq) {
774 <          if (iLoop == PAIR_LOOP) {
775 <            vij = 0.0;
776 <            fij.zero();
777 <            eField1.zero();
778 <            eField2.zero();
779 <          }
732 >          cg2 = neighborList_[m2];
733            
734 <          in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
735 <                                                     rgrp);
736 <
737 <          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
738 <          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
739 <
740 <          if (doHeatFlux_)
741 <            gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
742 <
743 <          for (ia = atomListRow.begin();
744 <               ia != atomListRow.end(); ++ia) {            
745 <            atom1 = (*ia);
746 <
747 <            for (jb = atomListColumn.begin();
748 <                 jb != atomListColumn.end(); ++jb) {              
749 <              atom2 = (*jb);
750 <
751 <              if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
752 <
753 <                vpair = 0.0;
754 <                workPot = 0.0;
755 <                exPot = 0.0;
756 <                f1.zero();
757 <                dVdFQ1 = 0.0;
758 <                dVdFQ2 = 0.0;
759 <
760 <                fDecomp_->fillInteractionData(idat, atom1, atom2);
761 <
762 <                topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
763 <                vdwMult = vdwScale_[topoDist];
764 <                electroMult = electrostaticScale_[topoDist];
765 <
766 <                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
767 <                  idat.d = &d_grp;
768 <                  idat.r2 = &rgrpsq;
769 <                  if (doHeatFlux_)
770 <                    vel2 = gvel2;
771 <                } else {
772 <                  d = fDecomp_->getInteratomicVector(atom1, atom2);
773 <                  curSnapshot->wrapVector( d );
774 <                  r2 = d.lengthSquare();
775 <                  idat.d = &d;
776 <                  idat.r2 = &r2;
777 <                  if (doHeatFlux_)
778 <                    vel2 = fDecomp_->getAtomVelocityColumn(atom2);
779 <                }
780 <              
781 <                r = sqrt( *(idat.r2) );
782 <                idat.rij = &r;
783 <
784 <                if (iLoop == PREPAIR_LOOP) {
785 <                  interactionMan_->doPrePair(idat);
786 <                } else {
787 <                  interactionMan_->doPair(idat);
788 <                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
789 <                  vij += vpair;
790 <                  fij += f1;
791 <                  stressTensor -= outProduct( *(idat.d), f1);
792 <                  if (doHeatFlux_)
793 <                    fDecomp_->addToHeatFlux(*(idat.d) * dot(f1, vel2));
734 >          d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
735 >        
736 >          // already wrapped in the getIntergroupVector call:
737 >          // curSnapshot->wrapVector(d_grp);        
738 >          rgrpsq = d_grp.lengthSquare();
739 >          
740 >          if (rgrpsq < rCutSq_) {
741 >            if (iLoop == PAIR_LOOP) {
742 >              vij = 0.0;
743 >              fij.zero();
744 >              eField1.zero();
745 >              eField2.zero();
746 >              sPot1 = 0.0;
747 >              sPot2 = 0.0;
748 >            }
749 >            
750 >            in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
751 >                                                       rgrp);
752 >            
753 >            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
754 >            
755 >            if (doHeatFlux_)
756 >              gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
757 >            
758 >            for (ia = atomListRow.begin();
759 >                 ia != atomListRow.end(); ++ia) {            
760 >              atom1 = (*ia);
761 >              
762 >              for (jb = atomListColumn.begin();
763 >                   jb != atomListColumn.end(); ++jb) {              
764 >                atom2 = (*jb);
765 >                
766 >                if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
767 >                  
768 >                  vpair = 0.0;
769 >                  workPot = 0.0;
770 >                  exPot = 0.0;
771 >                  f1.zero();
772 >                  dVdFQ1 = 0.0;
773 >                  dVdFQ2 = 0.0;
774 >                  
775 >                  fDecomp_->fillInteractionData(idat, atom1, atom2, newAtom1);
776 >                  
777 >                  topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
778 >                  vdwMult = vdwScale_[topoDist];
779 >                  electroMult = electrostaticScale_[topoDist];
780 >                  
781 >                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
782 >                    idat.d = &d_grp;
783 >                    idat.r2 = &rgrpsq;
784 >                    if (doHeatFlux_)
785 >                      vel2 = gvel2;
786 >                  } else {
787 >                    d = fDecomp_->getInteratomicVector(atom1, atom2);
788 >                    curSnapshot->wrapVector( d );
789 >                    r2 = d.lengthSquare();
790 >                    idat.d = &d;
791 >                    idat.r2 = &r2;
792 >                    if (doHeatFlux_)
793 >                      vel2 = fDecomp_->getAtomVelocityColumn(atom2);
794 >                  }
795 >                  
796 >                  r = sqrt( *(idat.r2) );
797 >                  idat.rij = &r;
798 >                  
799 >                  if (iLoop == PREPAIR_LOOP) {
800 >                    interactionMan_->doPrePair(idat);
801 >                  } else {
802 >                    interactionMan_->doPair(idat);
803 >                    fDecomp_->unpackInteractionData(idat, atom1, atom2);
804 >                    vij += vpair;
805 >                    fij += f1;
806 >                    stressTensor -= outProduct( *(idat.d), f1);
807 >                    if (doHeatFlux_)
808 >                      fDecomp_->addToHeatFlux(*(idat.d) * dot(f1, vel2));
809 >                  }
810                  }
811                }
812              }
813 <          }
814 <
815 <          if (iLoop == PAIR_LOOP) {
816 <            if (in_switching_region) {
817 <              swderiv = vij * dswdr / rgrp;
818 <              fg = swderiv * d_grp;
819 <              fij += fg;
820 <
821 <              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
822 <                if (!fDecomp_->skipAtomPair(atomListRow[0],
823 <                                            atomListColumn[0],
855 <                                            cg1, cg2)) {
813 >            
814 >            if (iLoop == PAIR_LOOP) {
815 >              if (in_switching_region) {
816 >                swderiv = vij * dswdr / rgrp;
817 >                fg = swderiv * d_grp;
818 >                fij += fg;
819 >                
820 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
821 >                  if (!fDecomp_->skipAtomPair(atomListRow[0],
822 >                                              atomListColumn[0],
823 >                                              cg1, cg2)) {
824                    stressTensor -= outProduct( *(idat.d), fg);
825                    if (doHeatFlux_)
826                      fDecomp_->addToHeatFlux(*(idat.d) * dot(fg, vel2));
827 <                }                
828 <              }
829 <          
830 <              for (ia = atomListRow.begin();
831 <                   ia != atomListRow.end(); ++ia) {            
832 <                atom1 = (*ia);                
833 <                mf = fDecomp_->getMassFactorRow(atom1);
834 <                // fg is the force on atom ia due to cutoff group's
835 <                // presence in switching region
836 <                fg = swderiv * d_grp * mf;
837 <                fDecomp_->addForceToAtomRow(atom1, fg);
838 <                if (atomListRow.size() > 1) {
839 <                  if (info_->usesAtomicVirial()) {
840 <                    // find the distance between the atom
841 <                    // and the center of the cutoff group:
842 <                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
843 <                    stressTensor -= outProduct(dag, fg);
844 <                    if (doHeatFlux_)
845 <                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
827 >                  }                
828 >                }
829 >                
830 >                for (ia = atomListRow.begin();
831 >                     ia != atomListRow.end(); ++ia) {            
832 >                  atom1 = (*ia);                
833 >                  mf = fDecomp_->getMassFactorRow(atom1);
834 >                  // fg is the force on atom ia due to cutoff group's
835 >                  // presence in switching region
836 >                  fg = swderiv * d_grp * mf;
837 >                  fDecomp_->addForceToAtomRow(atom1, fg);
838 >                  if (atomListRow.size() > 1) {
839 >                    if (info_->usesAtomicVirial()) {
840 >                      // find the distance between the atom
841 >                      // and the center of the cutoff group:
842 >                      dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
843 >                      stressTensor -= outProduct(dag, fg);
844 >                      if (doHeatFlux_)
845 >                        fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
846 >                    }
847                    }
848                  }
849 <              }
850 <              for (jb = atomListColumn.begin();
851 <                   jb != atomListColumn.end(); ++jb) {              
852 <                atom2 = (*jb);
853 <                mf = fDecomp_->getMassFactorColumn(atom2);
854 <                // fg is the force on atom jb due to cutoff group's
855 <                // presence in switching region
856 <                fg = -swderiv * d_grp * mf;
857 <                fDecomp_->addForceToAtomColumn(atom2, fg);
858 <
859 <                if (atomListColumn.size() > 1) {
860 <                  if (info_->usesAtomicVirial()) {
861 <                    // find the distance between the atom
862 <                    // and the center of the cutoff group:
863 <                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
864 <                    stressTensor -= outProduct(dag, fg);
865 <                    if (doHeatFlux_)
866 <                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
849 >                for (jb = atomListColumn.begin();
850 >                     jb != atomListColumn.end(); ++jb) {              
851 >                  atom2 = (*jb);
852 >                  mf = fDecomp_->getMassFactorColumn(atom2);
853 >                  // fg is the force on atom jb due to cutoff group's
854 >                  // presence in switching region
855 >                  fg = -swderiv * d_grp * mf;
856 >                  fDecomp_->addForceToAtomColumn(atom2, fg);
857 >                  
858 >                  if (atomListColumn.size() > 1) {
859 >                    if (info_->usesAtomicVirial()) {
860 >                      // find the distance between the atom
861 >                      // and the center of the cutoff group:
862 >                      dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
863 >                      stressTensor -= outProduct(dag, fg);
864 >                      if (doHeatFlux_)
865 >                        fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
866 >                    }
867                    }
868                  }
869                }
870 +              //if (!info_->usesAtomicVirial()) {
871 +              //  stressTensor -= outProduct(d_grp, fij);
872 +              //  if (doHeatFlux_)
873 +              //     fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
874 +              //}
875              }
902 –            //if (!info_->usesAtomicVirial()) {
903 –            //  stressTensor -= outProduct(d_grp, fij);
904 –            //  if (doHeatFlux_)
905 –            //     fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
906 –            //}
876            }
877          }
878 +        newAtom1 = false;
879        }
880 <
880 >        
881        if (iLoop == PREPAIR_LOOP) {
882          if (info_->requiresPrepair()) {
883 <
883 >          
884            fDecomp_->collectIntermediateData();
885 <
885 >          
886            for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
887              fDecomp_->fillSelfData(sdat, atom1);
888              interactionMan_->doPreForce(sdat);
889            }
890 <
890 >          
891            fDecomp_->distributeIntermediateData();
892 <
892 >          
893          }
894        }
895      }
# Line 948 | Line 918 | namespace OpenMD {
918      curSnapshot->setLongRangePotential(longRangePotential);
919      
920      curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) +
921 <                                         *(fDecomp_->getExcludedPotential()));
921 >                                       *(fDecomp_->getExcludedPotential()));
922  
923    }
924  
# Line 977 | Line 947 | namespace OpenMD {
947      }
948      
949   #ifdef IS_MPI
950 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9,
951 <                              MPI::REALTYPE, MPI::SUM);
950 >    MPI_Allreduce(MPI_IN_PLACE, stressTensor.getArrayPointer(), 9,
951 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
952   #endif
953      curSnapshot->setStressTensor(stressTensor);
954      
955      if (info_->getSimParams()->getUseLongRangeCorrections()) {
956        /*
957 <      RealType vol = curSnapshot->getVolume();
958 <      RealType Elrc(0.0);
959 <      RealType Wlrc(0.0);
957 >        RealType vol = curSnapshot->getVolume();
958 >        RealType Elrc(0.0);
959 >        RealType Wlrc(0.0);
960  
961 <      set<AtomType*>::iterator i;
962 <      set<AtomType*>::iterator j;
961 >        set<AtomType*>::iterator i;
962 >        set<AtomType*>::iterator j;
963      
964 <      RealType n_i, n_j;
965 <      RealType rho_i, rho_j;
966 <      pair<RealType, RealType> LRI;
964 >        RealType n_i, n_j;
965 >        RealType rho_i, rho_j;
966 >        pair<RealType, RealType> LRI;
967        
968 <      for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
968 >        for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
969          n_i = RealType(info_->getGlobalCountOfType(*i));
970          rho_i = n_i /  vol;
971          for (j = atomTypes_.begin(); j != atomTypes_.end(); ++j) {
972 <          n_j = RealType(info_->getGlobalCountOfType(*j));
973 <          rho_j = n_j / vol;
972 >        n_j = RealType(info_->getGlobalCountOfType(*j));
973 >        rho_j = n_j / vol;
974            
975 <          LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) );
975 >        LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) );
976  
977 <          Elrc += n_i   * rho_j * LRI.first;
978 <          Wlrc -= rho_i * rho_j * LRI.second;
977 >        Elrc += n_i   * rho_j * LRI.first;
978 >        Wlrc -= rho_i * rho_j * LRI.second;
979          }
980 <      }
981 <      Elrc *= 2.0 * NumericConstant::PI;
982 <      Wlrc *= 2.0 * NumericConstant::PI;
980 >        }
981 >        Elrc *= 2.0 * NumericConstant::PI;
982 >        Wlrc *= 2.0 * NumericConstant::PI;
983  
984 <      RealType lrp = curSnapshot->getLongRangePotential();
985 <      curSnapshot->setLongRangePotential(lrp + Elrc);
986 <      stressTensor += Wlrc * SquareMatrix3<RealType>::identity();
987 <      curSnapshot->setStressTensor(stressTensor);
984 >        RealType lrp = curSnapshot->getLongRangePotential();
985 >        curSnapshot->setLongRangePotential(lrp + Elrc);
986 >        stressTensor += Wlrc * SquareMatrix3<RealType>::identity();
987 >        curSnapshot->setStressTensor(stressTensor);
988        */
989      
990      }

Diff Legend

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