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

Diff Legend

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