--- trunk/src/brains/ForceManager.cpp	2013/06/16 15:15:42	1879
+++ trunk/src/brains/ForceManager.cpp	2014/04/29 17:32:31	1993
@@ -99,7 +99,7 @@ namespace OpenMD {
    *      Use the maximum suggested value that was found.
    *
    * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, TAYLOR_SHIFTED, 
-   *                        or SHIFTED_POTENTIAL)
+   *                        SHIFTED_POTENTIAL, or EWALD_FULL)
    *      If cutoffMethod was explicitly set, use that choice.
    *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
    *
@@ -172,6 +172,7 @@ namespace OpenMD {
     stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;    
     stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
     stringToCutoffMethod["TAYLOR_SHIFTED"] = TAYLOR_SHIFTED;
+    stringToCutoffMethod["EWALD_FULL"] = EWALD_FULL;
   
     if (simParams_->haveCutoffMethod()) {
       string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
@@ -182,7 +183,7 @@ namespace OpenMD {
                 "ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
                 "\tShould be one of: "
                 "HARD, SWITCHED, SHIFTED_POTENTIAL, TAYLOR_SHIFTED,\n"
-                "\tor SHIFTED_FORCE\n",
+                "\tSHIFTED_FORCE, or EWALD_FULL\n",
                 cutMeth.c_str());
         painCave.isFatal = 1;
         painCave.severity = OPENMD_ERROR;
@@ -228,13 +229,15 @@ namespace OpenMD {
             cutoffMethod_ = SHIFTED_FORCE;
           } else if (myMethod == "TAYLOR_SHIFTED") {
             cutoffMethod_ = TAYLOR_SHIFTED;
+          } else if (myMethod == "EWALD_FULL") {
+            cutoffMethod_ = EWALD_FULL;
           }
         
           if (simParams_->haveSwitchingRadius()) 
             rSwitch_ = simParams_->getSwitchingRadius();
 
           if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE" ||
-              myMethod == "TAYLOR_SHIFTED") {
+              myMethod == "TAYLOR_SHIFTED" || myMethod == "EWALD_FULL") {
             if (simParams_->haveSwitchingRadius()){
               sprintf(painCave.errMsg,
                       "ForceManager::setupCutoffs : DEPRECATED ERROR MESSAGE\n"
@@ -414,6 +417,7 @@ namespace OpenMD {
       if (doHeatFlux_) doParticlePot_ = true;
 
       doElectricField_ = info_->getSimParams()->getOutputElectricField();
+      doSitePotential_ = info_->getSimParams()->getOutputSitePotential();
    
     }
 
@@ -634,14 +638,15 @@ namespace OpenMD {
     // Collect from all nodes.  This should eventually be moved into a
     // SystemDecomposition, but this is a better place than in
     // Thermo to do the collection.
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE, 
-                              MPI::SUM);
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE, 
-                              MPI::SUM);
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1, 
-                              MPI::REALTYPE, MPI::SUM);
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1, 
-                              MPI::REALTYPE, MPI::SUM);
+
+    MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE, 
+                  MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE, 
+                  MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1, 
+                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1, 
+                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
 #endif
 
     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
@@ -663,6 +668,7 @@ namespace OpenMD {
     DataStorage* config = &(curSnapshot->atomData);
     DataStorage* cgConfig = &(curSnapshot->cgData);
 
+
     //calculate the center of mass of cutoff group
 
     SimInfo::MoleculeIterator mi;
@@ -695,7 +701,7 @@ namespace OpenMD {
     RealType vij;
     Vector3d fij, fg, f1;
     tuple3<RealType, RealType, RealType> cuts;
-    RealType rCutSq;
+    RealType rCut, rCutSq, rListSq;
     bool in_switching_region;
     RealType sw, dswdr, swderiv;
     vector<int> atomListColumn, atomListRow;
@@ -705,15 +711,20 @@ namespace OpenMD {
     RealType vpair;
     RealType dVdFQ1(0.0);
     RealType dVdFQ2(0.0);
-    Vector3d eField1(0.0);
-    Vector3d eField2(0.0);
     potVec longRangePotential(0.0);
+    RealType reciprocalPotential(0.0);
     potVec workPot(0.0);
     potVec exPot(0.0);
+    Vector3d eField1(0.0);
+    Vector3d eField2(0.0);
+    RealType sPot1(0.0);
+    RealType sPot2(0.0);
+                   
     vector<int>::iterator ia, jb;
 
     int loopStart, loopEnd;
-
+    
+    idat.rcut = &rCut;
     idat.vdwMult = &vdwMult;
     idat.electroMult = &electroMult;
     idat.pot = &workPot;
@@ -724,13 +735,16 @@ namespace OpenMD {
     idat.dVdFQ1 = &dVdFQ1;
     idat.dVdFQ2 = &dVdFQ2;
     idat.eField1 = &eField1;
-    idat.eField2 = &eField2;   
+    idat.eField2 = &eField2; 
+    idat.sPot1 = &sPot1;
+    idat.sPot2 = &sPot2;
     idat.f1 = &f1;
     idat.sw = &sw;
     idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
     idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE || cutoffMethod_ == TAYLOR_SHIFTED) ? true : false;
     idat.doParticlePot = doParticlePot_;
     idat.doElectricField = doElectricField_;
+    idat.doSitePotential = doSitePotential_;
     sdat.doParticlePot = doParticlePot_;
     
     loopEnd = PAIR_LOOP;
@@ -746,31 +760,32 @@ namespace OpenMD {
         if (update_nlist) {
           if (!usePeriodicBoundaryConditions_)
             Mat3x3d bbox = thermo->getBoundingBox();
-          neighborList = fDecomp_->buildNeighborList();
+          fDecomp_->buildNeighborList(neighborList_);
         }
       }
 
-      for (vector<pair<int, int> >::iterator it = neighborList.begin(); 
-             it != neighborList.end(); ++it) {
+      for (vector<pair<int, int> >::iterator it = neighborList_.begin(); 
+             it != neighborList_.end(); ++it) {
                 
         cg1 = (*it).first;
         cg2 = (*it).second;
         
-        cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
+        fDecomp_->getGroupCutoffs(cg1, cg2, rCut, rCutSq, rListSq);
 
         d_grp  = fDecomp_->getIntergroupVector(cg1, cg2);
 
-        curSnapshot->wrapVector(d_grp);        
+        // already wrapped in the getIntergroupVector call:
+        // curSnapshot->wrapVector(d_grp);        
         rgrpsq = d_grp.lengthSquare();
-        rCutSq = cuts.second;
 
         if (rgrpsq < rCutSq) {
-          idat.rcut = &cuts.first;
           if (iLoop == PAIR_LOOP) {
             vij = 0.0;
             fij.zero();
             eField1.zero();
             eField2.zero();
+            sPot1 = 0.0;
+            sPot2 = 0.0;
           }
           
           in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr, 
@@ -822,7 +837,7 @@ namespace OpenMD {
                
                 r = sqrt( *(idat.r2) );
                 idat.rij = &r;
-               
+
                 if (iLoop == PREPAIR_LOOP) {
                   interactionMan_->doPrePair(idat);
                 } else {
@@ -918,9 +933,14 @@ namespace OpenMD {
         }
       }
     }
-   
+    
     // collects pairwise information
     fDecomp_->collectData();
+    if (cutoffMethod_ == EWALD_FULL) {
+      interactionMan_->doReciprocalSpaceSum(reciprocalPotential);
+
+      curSnapshot->setReciprocalPotential(reciprocalPotential);
+    }
         
     if (info_->requiresSelfCorrection()) {
       for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) { 
@@ -942,7 +962,6 @@ namespace OpenMD {
 
   }
 
-  
   void ForceManager::postCalculation() {
 
     vector<Perturbation*>::iterator pi;
@@ -968,8 +987,8 @@ namespace OpenMD {
     }
     
 #ifdef IS_MPI
-    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9, 
-                              MPI::REALTYPE, MPI::SUM);
+    MPI_Allreduce(MPI_IN_PLACE, stressTensor.getArrayPointer(), 9, 
+                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
 #endif
     curSnapshot->setStressTensor(stressTensor);