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 1788 by gezelter, Wed Aug 29 20:17:07 2012 UTC vs.
Revision 2067 by gezelter, Thu Mar 5 15:35:37 2015 UTC

# Line 35 | Line 35
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).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 44 | Line 44
44   * @file ForceManager.cpp
45   * @author tlin
46   * @date 11/09/2004
47 – * @time 10:39am
47   * @version 1.0
48   */
49  
# Line 58 | 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 68 | Line 68 | namespace OpenMD {
68   using namespace std;
69   namespace OpenMD {
70    
71 <  ForceManager::ForceManager(SimInfo * info) : info_(info) {
71 >  ForceManager::ForceManager(SimInfo * info) : initialized_(false), info_(info),
72 >                                               switcher_(NULL) {
73      forceField_ = info_->getForceField();
74      interactionMan_ = new InteractionManager();
75      fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_);
76 +    thermo = new Thermo(info_);
77    }
78  
79 +  ForceManager::~ForceManager() {
80 +    perturbations_.clear();
81 +    
82 +    delete switcher_;
83 +    delete interactionMan_;
84 +    delete fDecomp_;
85 +    delete thermo;
86 +  }
87 +  
88    /**
89     * setupCutoffs
90     *
91 <   * Sets the values of cutoffRadius, switchingRadius, cutoffMethod,
81 <   * 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 88 | Line 98 | namespace OpenMD {
98     *      simulation for suggested cutoff values (e.g. 2.5 * sigma).
99     *      Use the maximum suggested value that was found.
100     *
101 <   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE,
102 <   *                        or SHIFTED_POTENTIAL)
101 >   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, TAYLOR_SHIFTED,
102 >   *                        SHIFTED_POTENTIAL, or EWALD_FULL)
103     *      If cutoffMethod was explicitly set, use that choice.
104     *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
105     *
96 –   * cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
97 –   *      If cutoffPolicy was explicitly set, use that choice.
98 –   *      If cutoffPolicy was not explicitly set, use TRADITIONAL
99 –   *
106     * switchingRadius : realType
107     *  If the cutoffMethod was set to SWITCHED:
108     *      If the switchingRadius was explicitly set, use that value
# Line 109 | Line 115 | namespace OpenMD {
115    void ForceManager::setupCutoffs() {
116      
117      Globals* simParams_ = info_->getSimParams();
112 –    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
118      int mdFileVersion;
119      rCut_ = 0.0; //Needs a value for a later max() call;  
120      
# Line 118 | Line 123 | namespace OpenMD {
123      else
124        mdFileVersion = 0;
125    
126 +    // We need the list of simulated atom types to figure out cutoffs
127 +    // as well as long range corrections.
128 +
129 +    set<AtomType*>::iterator i;
130 +    set<AtomType*> atomTypes_;
131 +    atomTypes_ = info_->getSimulatedAtomTypes();
132 +
133      if (simParams_->haveCutoffRadius()) {
134        rCut_ = simParams_->getCutoffRadius();
135      } else {      
# Line 132 | Line 144 | namespace OpenMD {
144          rCut_ = 12.0;
145        } else {
146          RealType thisCut;
147 <        set<AtomType*>::iterator i;
136 <        set<AtomType*> atomTypes;
137 <        atomTypes = info_->getSimulatedAtomTypes();        
138 <        for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
147 >        for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
148            thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
149            rCut_ = max(thisCut, rCut_);
150          }
# Line 149 | Line 158 | namespace OpenMD {
158        }
159      }
160  
161 <    fDecomp_->setUserCutoff(rCut_);
161 >    fDecomp_->setCutoffRadius(rCut_);
162      interactionMan_->setCutoffRadius(rCut_);
163 +    rCutSq_ = rCut_ * rCut_;
164  
165      map<string, CutoffMethod> stringToCutoffMethod;
166      stringToCutoffMethod["HARD"] = HARD;
167      stringToCutoffMethod["SWITCHED"] = SWITCHED;
168      stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;    
169      stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
170 +    stringToCutoffMethod["TAYLOR_SHIFTED"] = TAYLOR_SHIFTED;
171 +    stringToCutoffMethod["EWALD_FULL"] = EWALD_FULL;
172    
173      if (simParams_->haveCutoffMethod()) {
174        string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
# Line 166 | Line 178 | namespace OpenMD {
178          sprintf(painCave.errMsg,
179                  "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
180                  "\tShould be one of: "
181 <                "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
181 >                "HARD, SWITCHED, SHIFTED_POTENTIAL, TAYLOR_SHIFTED,\n"
182 >                "\tSHIFTED_FORCE, or EWALD_FULL\n",
183                  cutMeth.c_str());
184          painCave.isFatal = 1;
185          painCave.severity = OPENMD_ERROR;
# Line 210 | Line 223 | namespace OpenMD {
223              cutoffMethod_ = SHIFTED_POTENTIAL;
224            } else if (myMethod == "SHIFTED_FORCE") {
225              cutoffMethod_ = SHIFTED_FORCE;
226 +          } else if (myMethod == "TAYLOR_SHIFTED") {
227 +            cutoffMethod_ = TAYLOR_SHIFTED;
228 +          } else if (myMethod == "EWALD_FULL") {
229 +            cutoffMethod_ = EWALD_FULL;
230            }
231          
232            if (simParams_->haveSwitchingRadius())
233              rSwitch_ = simParams_->getSwitchingRadius();
234  
235 <          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
235 >          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE" ||
236 >              myMethod == "TAYLOR_SHIFTED" || myMethod == "EWALD_FULL") {
237              if (simParams_->haveSwitchingRadius()){
238                sprintf(painCave.errMsg,
239                        "ForceManager::setupCutoffs : DEPRECATED ERROR MESSAGE\n"
# Line 252 | Line 270 | namespace OpenMD {
270          }
271        }
272      }
255 –
256 –    map<string, CutoffPolicy> stringToCutoffPolicy;
257 –    stringToCutoffPolicy["MIX"] = MIX;
258 –    stringToCutoffPolicy["MAX"] = MAX;
259 –    stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;    
260 –
261 –    string cutPolicy;
262 –    if (forceFieldOptions_.haveCutoffPolicy()){
263 –      cutPolicy = forceFieldOptions_.getCutoffPolicy();
264 –    }else if (simParams_->haveCutoffPolicy()) {
265 –      cutPolicy = simParams_->getCutoffPolicy();
266 –    }
267 –
268 –    if (!cutPolicy.empty()){
269 –      toUpper(cutPolicy);
270 –      map<string, CutoffPolicy>::iterator i;
271 –      i = stringToCutoffPolicy.find(cutPolicy);
272 –
273 –      if (i == stringToCutoffPolicy.end()) {
274 –        sprintf(painCave.errMsg,
275 –                "ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
276 –                "\tShould be one of: "
277 –                "MIX, MAX, or TRADITIONAL\n",
278 –                cutPolicy.c_str());
279 –        painCave.isFatal = 1;
280 –        painCave.severity = OPENMD_ERROR;
281 –        simError();
282 –      } else {
283 –        cutoffPolicy_ = i->second;
284 –      }
285 –    } else {
286 –      sprintf(painCave.errMsg,
287 –              "ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
288 –              "\tOpenMD will use TRADITIONAL.\n");
289 –      painCave.isFatal = 0;
290 –      painCave.severity = OPENMD_INFO;
291 –      simError();
292 –      cutoffPolicy_ = TRADITIONAL;        
293 –    }
294 –
295 –    fDecomp_->setCutoffPolicy(cutoffPolicy_);
273          
274      // create the switching function object:
275  
# Line 370 | Line 347 | namespace OpenMD {
347      }
348      switcher_->setSwitchType(sft_);
349      switcher_->setSwitch(rSwitch_, rCut_);
373 –    interactionMan_->setSwitchingRadius(rSwitch_);
350    }
351  
376 –
377 –
378 –  
352    void ForceManager::initialize() {
353  
354      if (!info_->isTopologyDone()) {
# Line 384 | Line 357 | namespace OpenMD {
357        interactionMan_->setSimInfo(info_);
358        interactionMan_->initialize();
359  
360 <      // We want to delay the cutoffs until after the interaction
361 <      // manager has set up the atom-atom interactions so that we can
362 <      // query them for suggested cutoff values
360 >      //! We want to delay the cutoffs until after the interaction
361 >      //! manager has set up the atom-atom interactions so that we can
362 >      //! query them for suggested cutoff values
363        setupCutoffs();
364  
365        info_->prepareTopology();      
# Line 394 | Line 367 | namespace OpenMD {
367        doParticlePot_ = info_->getSimParams()->getOutputParticlePotential();
368        doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux();
369        if (doHeatFlux_) doParticlePot_ = true;
370 +
371 +      doElectricField_ = info_->getSimParams()->getOutputElectricField();
372 +      doSitePotential_ = info_->getSimParams()->getOutputSitePotential();
373    
374      }
375  
376      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
377      
378 <    // Force fields can set options on how to scale van der Waals and
379 <    // electrostatic interactions for atoms connected via bonds, bends
380 <    // and torsions in this case the topological distance between
381 <    // atoms is:
382 <    // 0 = topologically unconnected
383 <    // 1 = bonded together
384 <    // 2 = connected via a bend
385 <    // 3 = connected via a torsion
378 >    //! Force fields can set options on how to scale van der Waals and
379 >    //! electrostatic interactions for atoms connected via bonds, bends
380 >    //! and torsions in this case the topological distance between
381 >    //! atoms is:
382 >    //! 0 = topologically unconnected
383 >    //! 1 = bonded together
384 >    //! 2 = connected via a bend
385 >    //! 3 = connected via a torsion
386      
387      vdwScale_.reserve(4);
388      fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
# Line 424 | Line 400 | namespace OpenMD {
400      electrostaticScale_[2] = fopts.getelectrostatic13scale();
401      electrostaticScale_[3] = fopts.getelectrostatic14scale();    
402      
403 <    if (info_->getSimParams()->haveElectricField()) {
404 <      ElectricField* eField = new ElectricField(info_);
403 >    if (info_->getSimParams()->haveUniformField()) {
404 >      UniformField* eField = new UniformField(info_);
405        perturbations_.push_back(eField);
406      }
407 <
408 <    fDecomp_->distributeInitialData();
409 <
410 <    initialized_ = true;
411 <
407 >    if (info_->getSimParams()->haveUniformGradientStrength() ||
408 >        info_->getSimParams()->haveUniformGradientDirection1() ||
409 >        info_->getSimParams()->haveUniformGradientDirection2() ) {
410 >      UniformGradient* eGrad = new UniformGradient(info_);
411 >      perturbations_.push_back(eGrad);
412 >    }
413 >    
414 >    usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions();
415 >    
416 >    fDecomp_->distributeInitialData();
417 >    
418 >    initialized_ = true;
419 >    
420    }
421 <
421 >  
422    void ForceManager::calcForces() {
423      
424      if (!initialized_) initialize();
425 <
425 >    
426      preCalculation();  
427      shortRangeInteractions();
428      longRangeInteractions();
# Line 612 | Line 596 | namespace OpenMD {
596      // Collect from all nodes.  This should eventually be moved into a
597      // SystemDecomposition, but this is a better place than in
598      // Thermo to do the collection.
599 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE,
600 <                              MPI::SUM);
601 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE,
602 <                              MPI::SUM);
603 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1,
604 <                              MPI::REALTYPE, MPI::SUM);
605 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1,
606 <                              MPI::REALTYPE, MPI::SUM);
599 >
600 >    MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE,
601 >                  MPI_SUM, MPI_COMM_WORLD);
602 >    MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE,
603 >                  MPI_SUM, MPI_COMM_WORLD);
604 >    MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1,
605 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
606 >    MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1,
607 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
608   #endif
609  
610      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 637 | Line 622 | namespace OpenMD {
622    
623    void ForceManager::longRangeInteractions() {
624  
640 –
625      Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
626      DataStorage* config = &(curSnapshot->atomData);
627      DataStorage* cgConfig = &(curSnapshot->cgData);
# Line 648 | Line 632 | namespace OpenMD {
632      Molecule* mol;
633      Molecule::CutoffGroupIterator ci;
634      CutoffGroup* cg;
635 <
636 <    if(info_->getNCutoffGroups() > 0){      
635 >    
636 >    if(info_->getNCutoffGroups() != info_->getNAtoms()){
637        for (mol = info_->beginMolecule(mi); mol != NULL;
638             mol = info_->nextMolecule(mi)) {
639          for(cg = mol->beginCutoffGroup(ci); cg != NULL;
# Line 673 | Line 657 | namespace OpenMD {
657      RealType electroMult, vdwMult;
658      RealType vij;
659      Vector3d fij, fg, f1;
676 –    tuple3<RealType, RealType, RealType> cuts;
677 –    RealType rCutSq;
660      bool in_switching_region;
661      RealType sw, dswdr, swderiv;
662 <    vector<int> atomListColumn, atomListRow, atomListLocal;
662 >    vector<int> atomListColumn, atomListRow;
663      InteractionData idat;
664      SelfData sdat;
665      RealType mf;
666      RealType vpair;
667      RealType dVdFQ1(0.0);
668      RealType dVdFQ2(0.0);
687 –    Vector3d eField1(0.0);
688 –    Vector3d eField2(0.0);
669      potVec longRangePotential(0.0);
670 +    RealType reciprocalPotential(0.0);
671      potVec workPot(0.0);
672      potVec exPot(0.0);
673 +    Vector3d eField1(0.0);
674 +    Vector3d eField2(0.0);
675 +    RealType sPot1(0.0);
676 +    RealType sPot2(0.0);
677 +    bool newAtom1;
678 +                  
679      vector<int>::iterator ia, jb;
680  
681      int loopStart, loopEnd;
682 <
682 >    
683 >    idat.rcut = &rCut_;
684      idat.vdwMult = &vdwMult;
685      idat.electroMult = &electroMult;
686      idat.pot = &workPot;
# Line 703 | Line 691 | namespace OpenMD {
691      idat.dVdFQ1 = &dVdFQ1;
692      idat.dVdFQ2 = &dVdFQ2;
693      idat.eField1 = &eField1;
694 <    idat.eField2 = &eField2;
694 >    idat.eField2 = &eField2;
695 >    idat.sPot1 = &sPot1;
696 >    idat.sPot2 = &sPot2;
697      idat.f1 = &f1;
698      idat.sw = &sw;
699      idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
700 <    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
700 >    idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE ||
701 >                         cutoffMethod_ == TAYLOR_SHIFTED) ? true : false;
702      idat.doParticlePot = doParticlePot_;
703 +    idat.doElectricField = doElectricField_;
704 +    idat.doSitePotential = doSitePotential_;
705      sdat.doParticlePot = doParticlePot_;
706      
707      loopEnd = PAIR_LOOP;
# Line 721 | Line 714 | namespace OpenMD {
714      
715        if (iLoop == loopStart) {
716          bool update_nlist = fDecomp_->checkNeighborList();
717 <        if (update_nlist)
718 <          neighborList = fDecomp_->buildNeighborList();
719 <      }            
720 <
721 <      for (vector<pair<int, int> >::iterator it = neighborList.begin();
722 <             it != neighborList.end(); ++it) {
730 <                
731 <        cg1 = (*it).first;
732 <        cg2 = (*it).second;
733 <        
734 <        cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
717 >        if (update_nlist) {
718 >          if (!usePeriodicBoundaryConditions_)
719 >            Mat3x3d bbox = thermo->getBoundingBox();
720 >          fDecomp_->buildNeighborList(neighborList_, point_);
721 >        }
722 >      }
723  
724 <        d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
724 >      for (cg1 = 0; cg1 < int(point_.size()) - 1; cg1++) {
725 >        
726 >        atomListRow = fDecomp_->getAtomsInGroupRow(cg1);        
727 >        newAtom1 = true;
728 >        
729 >        for (int m2 = point_[cg1]; m2 < point_[cg1+1]; m2++) {
730  
731 <        curSnapshot->wrapVector(d_grp);        
739 <        rgrpsq = d_grp.lengthSquare();
740 <        rCutSq = cuts.second;
741 <
742 <        if (rgrpsq < rCutSq) {
743 <          idat.rcut = &cuts.first;
744 <          if (iLoop == PAIR_LOOP) {
745 <            vij = 0.0;
746 <            fij = V3Zero;
747 <          }
731 >          cg2 = neighborList_[m2];
732            
733 <          in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
734 <                                                     rgrp);
735 <
736 <          atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
737 <          atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
738 <
739 <          if (doHeatFlux_)
740 <            gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
741 <
742 <          for (ia = atomListRow.begin();
743 <               ia != atomListRow.end(); ++ia) {            
744 <            atom1 = (*ia);
745 <
746 <            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 = V3Zero;
772 <                dVdFQ1 = 0.0;
773 <                dVdFQ2 = 0.0;
774 <
775 <                fDecomp_->fillInteractionData(idat, atom1, atom2);
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 <              }
733 >          d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
734 >        
735 >          // already wrapped in the getIntergroupVector call:
736 >          // curSnapshot->wrapVector(d_grp);        
737 >          rgrpsq = d_grp.lengthSquare();
738 >          
739 >          if (rgrpsq < rCutSq_) {
740 >            if (iLoop == PAIR_LOOP) {
741 >              vij = 0.0;
742 >              fij.zero();
743 >              eField1.zero();
744 >              eField2.zero();
745 >              sPot1 = 0.0;
746 >              sPot2 = 0.0;
747              }
748 <          }
749 <
750 <          if (iLoop == PAIR_LOOP) {
751 <            if (in_switching_region) {
752 <              swderiv = vij * dswdr / rgrp;
753 <              fg = swderiv * d_grp;
754 <              fij += fg;
755 <
756 <              if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
757 <                stressTensor -= outProduct( *(idat.d), fg);
758 <                if (doHeatFlux_)
759 <                  fDecomp_->addToHeatFlux(*(idat.d) * dot(fg, vel2));
748 >            
749 >            in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
750 >                                                       rgrp);
751 >            
752 >            atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
753 >            
754 >            if (doHeatFlux_)
755 >              gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
756 >            
757 >            for (ia = atomListRow.begin();
758 >                 ia != atomListRow.end(); ++ia) {            
759 >              atom1 = (*ia);
760 >              
761 >              for (jb = atomListColumn.begin();
762 >                   jb != atomListColumn.end(); ++jb) {              
763 >                atom2 = (*jb);
764                  
765 <              }
766 <          
767 <              for (ia = atomListRow.begin();
768 <                   ia != atomListRow.end(); ++ia) {            
769 <                atom1 = (*ia);                
770 <                mf = fDecomp_->getMassFactorRow(atom1);
771 <                // fg is the force on atom ia due to cutoff group's
772 <                // presence in switching region
773 <                fg = swderiv * d_grp * mf;
774 <                fDecomp_->addForceToAtomRow(atom1, fg);
775 <                if (atomListRow.size() > 1) {
776 <                  if (info_->usesAtomicVirial()) {
777 <                    // find the distance between the atom
778 <                    // and the center of the cutoff group:
779 <                    dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
780 <                    stressTensor -= outProduct(dag, fg);
765 >                if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
766 >                  
767 >                  vpair = 0.0;
768 >                  workPot = 0.0;
769 >                  exPot = 0.0;
770 >                  f1.zero();
771 >                  dVdFQ1 = 0.0;
772 >                  dVdFQ2 = 0.0;
773 >                  
774 >                  fDecomp_->fillInteractionData(idat, atom1, atom2, newAtom1);
775 >                  
776 >                  topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
777 >                  vdwMult = vdwScale_[topoDist];
778 >                  electroMult = electrostaticScale_[topoDist];
779 >                  
780 >                  if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
781 >                    idat.d = &d_grp;
782 >                    idat.r2 = &rgrpsq;
783                      if (doHeatFlux_)
784 <                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
784 >                      vel2 = gvel2;
785 >                  } else {
786 >                    d = fDecomp_->getInteratomicVector(atom1, atom2);
787 >                    curSnapshot->wrapVector( d );
788 >                    r2 = d.lengthSquare();
789 >                    idat.d = &d;
790 >                    idat.r2 = &r2;
791 >                    if (doHeatFlux_)
792 >                      vel2 = fDecomp_->getAtomVelocityColumn(atom2);
793                    }
794 +                  
795 +                  r = sqrt( *(idat.r2) );
796 +                  idat.rij = &r;
797 +                  
798 +                  if (iLoop == PREPAIR_LOOP) {
799 +                    interactionMan_->doPrePair(idat);
800 +                  } else {
801 +                    interactionMan_->doPair(idat);
802 +                    fDecomp_->unpackInteractionData(idat, atom1, atom2);
803 +                    vij += vpair;
804 +                    fij += f1;
805 +                    stressTensor -= outProduct( *(idat.d), f1);
806 +                    if (doHeatFlux_)
807 +                      fDecomp_->addToHeatFlux(*(idat.d) * dot(f1, vel2));
808 +                  }
809                  }
810                }
811 <              for (jb = atomListColumn.begin();
812 <                   jb != atomListColumn.end(); ++jb) {              
813 <                atom2 = (*jb);
814 <                mf = fDecomp_->getMassFactorColumn(atom2);
815 <                // fg is the force on atom jb due to cutoff group's
816 <                // presence in switching region
817 <                fg = -swderiv * d_grp * mf;
818 <                fDecomp_->addForceToAtomColumn(atom2, fg);
819 <
820 <                if (atomListColumn.size() > 1) {
821 <                  if (info_->usesAtomicVirial()) {
822 <                    // find the distance between the atom
823 <                    // and the center of the cutoff group:
824 <                    dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
825 <                    stressTensor -= outProduct(dag, fg);
826 <                    if (doHeatFlux_)
827 <                      fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
811 >            }
812 >            
813 >            if (iLoop == PAIR_LOOP) {
814 >              if (in_switching_region) {
815 >                swderiv = vij * dswdr / rgrp;
816 >                fg = swderiv * d_grp;
817 >                fij += fg;
818 >                
819 >                if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
820 >                  if (!fDecomp_->skipAtomPair(atomListRow[0],
821 >                                              atomListColumn[0],
822 >                                              cg1, cg2)) {
823 >                  stressTensor -= outProduct( *(idat.d), fg);
824 >                  if (doHeatFlux_)
825 >                    fDecomp_->addToHeatFlux(*(idat.d) * dot(fg, vel2));
826 >                  }                
827 >                }
828 >                
829 >                for (ia = atomListRow.begin();
830 >                     ia != atomListRow.end(); ++ia) {            
831 >                  atom1 = (*ia);                
832 >                  mf = fDecomp_->getMassFactorRow(atom1);
833 >                  // fg is the force on atom ia due to cutoff group's
834 >                  // presence in switching region
835 >                  fg = swderiv * d_grp * mf;
836 >                  fDecomp_->addForceToAtomRow(atom1, fg);
837 >                  if (atomListRow.size() > 1) {
838 >                    if (info_->usesAtomicVirial()) {
839 >                      // find the distance between the atom
840 >                      // and the center of the cutoff group:
841 >                      dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
842 >                      stressTensor -= outProduct(dag, fg);
843 >                      if (doHeatFlux_)
844 >                        fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
845 >                    }
846                    }
847                  }
848 +                for (jb = atomListColumn.begin();
849 +                     jb != atomListColumn.end(); ++jb) {              
850 +                  atom2 = (*jb);
851 +                  mf = fDecomp_->getMassFactorColumn(atom2);
852 +                  // fg is the force on atom jb due to cutoff group's
853 +                  // presence in switching region
854 +                  fg = -swderiv * d_grp * mf;
855 +                  fDecomp_->addForceToAtomColumn(atom2, fg);
856 +                  
857 +                  if (atomListColumn.size() > 1) {
858 +                    if (info_->usesAtomicVirial()) {
859 +                      // find the distance between the atom
860 +                      // and the center of the cutoff group:
861 +                      dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
862 +                      stressTensor -= outProduct(dag, fg);
863 +                      if (doHeatFlux_)
864 +                        fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
865 +                    }
866 +                  }
867 +                }
868                }
869 +              //if (!info_->usesAtomicVirial()) {
870 +              //  stressTensor -= outProduct(d_grp, fij);
871 +              //  if (doHeatFlux_)
872 +              //     fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
873 +              //}
874              }
867 –            //if (!info_->usesAtomicVirial()) {
868 –            //  stressTensor -= outProduct(d_grp, fij);
869 –            //  if (doHeatFlux_)
870 –            //     fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
871 –            //}
875            }
876          }
877 +        newAtom1 = false;
878        }
879 <
879 >        
880        if (iLoop == PREPAIR_LOOP) {
881          if (info_->requiresPrepair()) {
882 <
882 >          
883            fDecomp_->collectIntermediateData();
884 <
884 >          
885            for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
886              fDecomp_->fillSelfData(sdat, atom1);
887              interactionMan_->doPreForce(sdat);
888            }
889 <
889 >          
890            fDecomp_->distributeIntermediateData();
891 <
891 >          
892          }
893        }
894      }
895 <  
895 >    
896      // collects pairwise information
897      fDecomp_->collectData();
898 +    if (cutoffMethod_ == EWALD_FULL) {
899 +      interactionMan_->doReciprocalSpaceSum(reciprocalPotential);
900 +
901 +      curSnapshot->setReciprocalPotential(reciprocalPotential);
902 +    }
903          
904      if (info_->requiresSelfCorrection()) {
905        for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
# Line 908 | Line 917 | namespace OpenMD {
917      curSnapshot->setLongRangePotential(longRangePotential);
918      
919      curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) +
920 <                                         *(fDecomp_->getExcludedPotential()));
920 >                                       *(fDecomp_->getExcludedPotential()));
921  
922    }
923  
915 –  
924    void ForceManager::postCalculation() {
925  
926      vector<Perturbation*>::iterator pi;
# Line 938 | Line 946 | namespace OpenMD {
946      }
947      
948   #ifdef IS_MPI
949 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9,
950 <                              MPI::REALTYPE, MPI::SUM);
949 >    MPI_Allreduce(MPI_IN_PLACE, stressTensor.getArrayPointer(), 9,
950 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
951   #endif
952      curSnapshot->setStressTensor(stressTensor);
953      
954 +    if (info_->getSimParams()->getUseLongRangeCorrections()) {
955 +      /*
956 +        RealType vol = curSnapshot->getVolume();
957 +        RealType Elrc(0.0);
958 +        RealType Wlrc(0.0);
959 +
960 +        set<AtomType*>::iterator i;
961 +        set<AtomType*>::iterator j;
962 +    
963 +        RealType n_i, n_j;
964 +        RealType rho_i, rho_j;
965 +        pair<RealType, RealType> LRI;
966 +      
967 +        for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
968 +        n_i = RealType(info_->getGlobalCountOfType(*i));
969 +        rho_i = n_i /  vol;
970 +        for (j = atomTypes_.begin(); j != atomTypes_.end(); ++j) {
971 +        n_j = RealType(info_->getGlobalCountOfType(*j));
972 +        rho_j = n_j / vol;
973 +          
974 +        LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) );
975 +
976 +        Elrc += n_i   * rho_j * LRI.first;
977 +        Wlrc -= rho_i * rho_j * LRI.second;
978 +        }
979 +        }
980 +        Elrc *= 2.0 * NumericConstant::PI;
981 +        Wlrc *= 2.0 * NumericConstant::PI;
982 +
983 +        RealType lrp = curSnapshot->getLongRangePotential();
984 +        curSnapshot->setLongRangePotential(lrp + Elrc);
985 +        stressTensor += Wlrc * SquareMatrix3<RealType>::identity();
986 +        curSnapshot->setStressTensor(stressTensor);
987 +      */
988 +    
989 +    }
990    }
991 < } //end namespace OpenMD
991 > }

Diff Legend

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