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 2047 by gezelter, Thu Dec 11 16:16:43 2014 UTC vs.
Revision 2057 by gezelter, Tue Mar 3 15:22:26 2015 UTC

# Line 88 | Line 88 | namespace OpenMD {
88    /**
89     * setupCutoffs
90     *
91 <   * Sets the values of cutoffRadius, switchingRadius, cutoffMethod,
92 <   * 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 104 | Line 103 | namespace OpenMD {
103     *      If cutoffMethod was explicitly set, use that choice.
104     *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
105     *
107 –   * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
108 –   *      If cutoffPolicy was explicitly set, use that choice.
109 –   *      If cutoffPolicy was not explicitly set, use TRADITIONAL
110 –   *
106     * switchingRadius : realType
107     *  If the cutoffMethod was set to SWITCHED:
108     *      If the switchingRadius was explicitly set, use that value
# Line 164 | 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 269 | namespace OpenMD {
269              }
270            }
271          }
276 –      }
277 –    }
278 –
279 –    map<string, CutoffPolicy> stringToCutoffPolicy;
280 –    stringToCutoffPolicy["MIX"] = MIX;
281 –    stringToCutoffPolicy["MAX"] = MAX;
282 –    stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;    
283 –
284 –    string cutPolicy;
285 –    if (forceFieldOptions_.haveCutoffPolicy()){
286 –      cutPolicy = forceFieldOptions_.getCutoffPolicy();
287 –    }else if (simParams_->haveCutoffPolicy()) {
288 –      cutPolicy = simParams_->getCutoffPolicy();
289 –    }
290 –
291 –    if (!cutPolicy.empty()){
292 –      toUpper(cutPolicy);
293 –      map<string, CutoffPolicy>::iterator i;
294 –      i = stringToCutoffPolicy.find(cutPolicy);
295 –
296 –      if (i == stringToCutoffPolicy.end()) {
297 –        sprintf(painCave.errMsg,
298 –                "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
299 –                "\tShould be one of: "
300 –                "MIX, MAX, or TRADITIONAL\n",
301 –                cutPolicy.c_str());
302 –        painCave.isFatal = 1;
303 –        painCave.severity = OPENMD_ERROR;
304 –        simError();
305 –      } else {
306 –        cutoffPolicy_ = i->second;
272        }
308 –    } else {
309 –      sprintf(painCave.errMsg,
310 –              "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
311 –              "\tOpenMD will use TRADITIONAL.\n");
312 –      painCave.isFatal = 0;
313 –      painCave.severity = OPENMD_INFO;
314 –      simError();
315 –      cutoffPolicy_ = TRADITIONAL;        
273      }
317 –
318 –    fDecomp_->setCutoffPolicy(cutoffPolicy_);
274          
275      // create the switching function object:
276  
# Line 671 | Line 626 | namespace OpenMD {
626      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
627      DataStorage* config = &(curSnapshot->atomData);
628      DataStorage* cgConfig = &(curSnapshot->cgData);
629 <
629 >    int jstart, jend;
630  
631      //calculate the center of mass of cutoff group
632  
# Line 679 | 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 723 | Line 678 | namespace OpenMD {
678      Vector3d eField2(0.0);
679      RealType sPot1(0.0);
680      RealType sPot2(0.0);
681 +    bool newAtom1;
682                    
683      vector<int>::iterator ia, jb;
684  
685      int loopStart, loopEnd;
686      
687 <    idat.rcut = &rCut;
687 >    idat.rcut = &rCut_;
688      idat.vdwMult = &vdwMult;
689      idat.electroMult = &electroMult;
690      idat.pot = &workPot;
# Line 765 | Line 721 | namespace OpenMD {
721          if (update_nlist) {
722            if (!usePeriodicBoundaryConditions_)
723              Mat3x3d bbox = thermo->getBoundingBox();
724 <          fDecomp_->buildNeighborList(neighborList_);
724 >          fDecomp_->buildNeighborList(neighborList_, point_);
725          }
726        }
727  
728 <      for (vector<pair<int, int> >::iterator it = neighborList_.begin();
773 <           it != neighborList_.end(); ++it) {
774 <                
775 <        cg1 = (*it).first;
776 <        cg2 = (*it).second;
728 >      for (unsigned int cg1 = 0; cg1 < point_.size() - 1; cg1++) {
729          
730 <        fDecomp_->getGroupCutoffs(cg1, cg2, rCut, rCutSq, rListSq);
731 <
732 <        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
733 <
782 <        // already wrapped in the getIntergroupVector call:
783 <        // curSnapshot->wrapVector(d_grp);        
784 <        rgrpsq = d_grp.lengthSquare();
785 <
786 <        if (rgrpsq < rCutSq) {
787 <          if (iLoop == PAIR_LOOP) {
788 <            vij = 0.0;
789 <            fij.zero();
790 <            eField1.zero();
791 <            eField2.zero();
792 <            sPot1 = 0.0;
793 <            sPot2 = 0.0;
794 <          }
795 <          
796 <          in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
797 <                                                     rgrp);
730 >        atomListRow = fDecomp_->getAtomsInGroupRow(cg1);        
731 >        newAtom1 = true;
732 >        
733 >        for (int m2 = point_[cg1]; m2 < point_[cg1+1]; m2++) {
734  
735 <          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
736 <          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
737 <
738 <          if (doHeatFlux_)
739 <            gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
740 <
741 <          for (ia = atomListRow.begin();
742 <               ia != atomListRow.end(); ++ia) {            
743 <            atom1 = (*ia);
744 <
745 <            for (jb = atomListColumn.begin();
746 <                 jb != atomListColumn.end(); ++jb) {              
747 <              atom2 = (*jb);
748 <
749 <              if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
750 <
751 <                vpair = 0.0;
752 <                workPot = 0.0;
753 <                exPot = 0.0;
754 <                f1.zero();
755 <                dVdFQ1 = 0.0;
756 <                dVdFQ2 = 0.0;
757 <
758 <                fDecomp_->fillInteractionData(idat, atom1, atom2);
759 <
760 <                topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
761 <                vdwMult = vdwScale_[topoDist];
762 <                electroMult = electrostaticScale_[topoDist];
763 <
764 <                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
765 <                  idat.d = &d_grp;
766 <                  idat.r2 = &rgrpsq;
767 <                  if (doHeatFlux_)
768 <                    vel2 = gvel2;
769 <                } else {
770 <                  d = fDecomp_->getInteratomicVector(atom1, atom2);
771 <                  curSnapshot->wrapVector( d );
772 <                  r2 = d.lengthSquare();
773 <                  idat.d = &d;
774 <                  idat.r2 = &r2;
775 <                  if (doHeatFlux_)
776 <                    vel2 = fDecomp_->getAtomVelocityColumn(atom2);
777 <                }
778 <              
779 <                r = sqrt( *(idat.r2) );
780 <                idat.rij = &r;
781 <
782 <                if (iLoop == PREPAIR_LOOP) {
783 <                  interactionMan_->doPrePair(idat);
784 <                } else {
785 <                  interactionMan_->doPair(idat);
786 <                  fDecomp_->unpackInteractionData(idat, atom1, atom2);
787 <                  vij += vpair;
788 <                  fij += f1;
789 <                  stressTensor -= outProduct( *(idat.d), f1);
790 <                  if (doHeatFlux_)
791 <                    fDecomp_->addToHeatFlux(*(idat.d) * dot(f1, vel2));
735 >          cg2 = neighborList_[m2];
736 >          
737 >          d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
738 >        
739 >          // already wrapped in the getIntergroupVector call:
740 >          // curSnapshot->wrapVector(d_grp);        
741 >          rgrpsq = d_grp.lengthSquare();
742 >          
743 >          if (rgrpsq < rCutSq_) {
744 >            if (iLoop == PAIR_LOOP) {
745 >              vij = 0.0;
746 >              fij.zero();
747 >              eField1.zero();
748 >              eField2.zero();
749 >              sPot1 = 0.0;
750 >              sPot2 = 0.0;
751 >            }
752 >            
753 >            in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
754 >                                                       rgrp);
755 >            
756 >            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
757 >            
758 >            if (doHeatFlux_)
759 >              gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
760 >            
761 >            for (ia = atomListRow.begin();
762 >                 ia != atomListRow.end(); ++ia) {            
763 >              atom1 = (*ia);
764 >              
765 >              for (jb = atomListColumn.begin();
766 >                   jb != atomListColumn.end(); ++jb) {              
767 >                atom2 = (*jb);
768 >                
769 >                if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
770 >                  
771 >                  vpair = 0.0;
772 >                  workPot = 0.0;
773 >                  exPot = 0.0;
774 >                  f1.zero();
775 >                  dVdFQ1 = 0.0;
776 >                  dVdFQ2 = 0.0;
777 >                  
778 >                  fDecomp_->fillInteractionData(idat, atom1, atom2, newAtom1);
779 >                  
780 >                  topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
781 >                  vdwMult = vdwScale_[topoDist];
782 >                  electroMult = electrostaticScale_[topoDist];
783 >                  
784 >                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
785 >                    idat.d = &d_grp;
786 >                    idat.r2 = &rgrpsq;
787 >                    if (doHeatFlux_)
788 >                      vel2 = gvel2;
789 >                  } else {
790 >                    d = fDecomp_->getInteratomicVector(atom1, atom2);
791 >                    curSnapshot->wrapVector( d );
792 >                    r2 = d.lengthSquare();
793 >                    idat.d = &d;
794 >                    idat.r2 = &r2;
795 >                    if (doHeatFlux_)
796 >                      vel2 = fDecomp_->getAtomVelocityColumn(atom2);
797 >                  }
798 >                  
799 >                  r = sqrt( *(idat.r2) );
800 >                  idat.rij = &r;
801 >                  
802 >                  if (iLoop == PREPAIR_LOOP) {
803 >                    interactionMan_->doPrePair(idat);
804 >                  } else {
805 >                    interactionMan_->doPair(idat);
806 >                    fDecomp_->unpackInteractionData(idat, atom1, atom2);
807 >                    vij += vpair;
808 >                    fij += f1;
809 >                    stressTensor -= outProduct( *(idat.d), f1);
810 >                    if (doHeatFlux_)
811 >                      fDecomp_->addToHeatFlux(*(idat.d) * dot(f1, vel2));
812 >                  }
813                  }
814                }
815              }
816 <          }
817 <
818 <          if (iLoop == PAIR_LOOP) {
819 <            if (in_switching_region) {
820 <              swderiv = vij * dswdr / rgrp;
821 <              fg = swderiv * d_grp;
822 <              fij += fg;
823 <
824 <              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
825 <                if (!fDecomp_->skipAtomPair(atomListRow[0],
826 <                                            atomListColumn[0],
870 <                                            cg1, cg2)) {
816 >            
817 >            if (iLoop == PAIR_LOOP) {
818 >              if (in_switching_region) {
819 >                swderiv = vij * dswdr / rgrp;
820 >                fg = swderiv * d_grp;
821 >                fij += fg;
822 >                
823 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
824 >                  if (!fDecomp_->skipAtomPair(atomListRow[0],
825 >                                              atomListColumn[0],
826 >                                              cg1, cg2)) {
827                    stressTensor -= outProduct( *(idat.d), fg);
828                    if (doHeatFlux_)
829                      fDecomp_->addToHeatFlux(*(idat.d) * dot(fg, vel2));
830 <                }                
831 <              }
832 <          
833 <              for (ia = atomListRow.begin();
834 <                   ia != atomListRow.end(); ++ia) {            
835 <                atom1 = (*ia);                
836 <                mf = fDecomp_->getMassFactorRow(atom1);
837 <                // fg is the force on atom ia due to cutoff group's
838 <                // presence in switching region
839 <                fg = swderiv * d_grp * mf;
840 <                fDecomp_->addForceToAtomRow(atom1, fg);
841 <                if (atomListRow.size() > 1) {
842 <                  if (info_->usesAtomicVirial()) {
843 <                    // find the distance between the atom
844 <                    // and the center of the cutoff group:
845 <                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
846 <                    stressTensor -= outProduct(dag, fg);
847 <                    if (doHeatFlux_)
848 <                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
830 >                  }                
831 >                }
832 >                
833 >                for (ia = atomListRow.begin();
834 >                     ia != atomListRow.end(); ++ia) {            
835 >                  atom1 = (*ia);                
836 >                  mf = fDecomp_->getMassFactorRow(atom1);
837 >                  // fg is the force on atom ia due to cutoff group's
838 >                  // presence in switching region
839 >                  fg = swderiv * d_grp * mf;
840 >                  fDecomp_->addForceToAtomRow(atom1, fg);
841 >                  if (atomListRow.size() > 1) {
842 >                    if (info_->usesAtomicVirial()) {
843 >                      // find the distance between the atom
844 >                      // and the center of the cutoff group:
845 >                      dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
846 >                      stressTensor -= outProduct(dag, fg);
847 >                      if (doHeatFlux_)
848 >                        fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
849 >                    }
850                    }
851                  }
852 <              }
853 <              for (jb = atomListColumn.begin();
854 <                   jb != atomListColumn.end(); ++jb) {              
855 <                atom2 = (*jb);
856 <                mf = fDecomp_->getMassFactorColumn(atom2);
857 <                // fg is the force on atom jb due to cutoff group's
858 <                // presence in switching region
859 <                fg = -swderiv * d_grp * mf;
860 <                fDecomp_->addForceToAtomColumn(atom2, fg);
861 <
862 <                if (atomListColumn.size() > 1) {
863 <                  if (info_->usesAtomicVirial()) {
864 <                    // find the distance between the atom
865 <                    // and the center of the cutoff group:
866 <                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
867 <                    stressTensor -= outProduct(dag, fg);
868 <                    if (doHeatFlux_)
869 <                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
852 >                for (jb = atomListColumn.begin();
853 >                     jb != atomListColumn.end(); ++jb) {              
854 >                  atom2 = (*jb);
855 >                  mf = fDecomp_->getMassFactorColumn(atom2);
856 >                  // fg is the force on atom jb due to cutoff group's
857 >                  // presence in switching region
858 >                  fg = -swderiv * d_grp * mf;
859 >                  fDecomp_->addForceToAtomColumn(atom2, fg);
860 >                  
861 >                  if (atomListColumn.size() > 1) {
862 >                    if (info_->usesAtomicVirial()) {
863 >                      // find the distance between the atom
864 >                      // and the center of the cutoff group:
865 >                      dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
866 >                      stressTensor -= outProduct(dag, fg);
867 >                      if (doHeatFlux_)
868 >                        fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
869 >                    }
870                    }
871                  }
872                }
873 +              //if (!info_->usesAtomicVirial()) {
874 +              //  stressTensor -= outProduct(d_grp, fij);
875 +              //  if (doHeatFlux_)
876 +              //     fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
877 +              //}
878              }
917 –            //if (!info_->usesAtomicVirial()) {
918 –            //  stressTensor -= outProduct(d_grp, fij);
919 –            //  if (doHeatFlux_)
920 –            //     fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
921 –            //}
879            }
880          }
881 +        newAtom1 = false;
882        }
883 <
883 >        
884        if (iLoop == PREPAIR_LOOP) {
885          if (info_->requiresPrepair()) {
886 <
886 >          
887            fDecomp_->collectIntermediateData();
888 <
888 >          
889            for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
890              fDecomp_->fillSelfData(sdat, atom1);
891              interactionMan_->doPreForce(sdat);
892            }
893 <
893 >          
894            fDecomp_->distributeIntermediateData();
895 <
895 >          
896          }
897        }
898      }

Diff Legend

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