Merge branch 'master' into alpha

plugins/lammps/USER-MBX/pair_mbx.cpp

@@ -200,6 +200,29 @@ void PairMBX::compute(int eflag, int vflag) {
 #endif
     }
 
+    if (mbx_parallel) {
+#ifdef _DEBUG
+        printf("[MBX] (%i) -- Computing electrostatics parallel\n", me);
+#endif
+
+        fix_mbx->mbxt_start(MBXT_ELE);
+        mbx_ele = ptr_mbx_local->ElectrostaticsMPIlocal(true, true);
+        fix_mbx->mbxt_stop(MBXT_ELE);
+        accumulate_f_local(true);
+
+    } else {
+#ifdef _DEBUG
+        printf("[MBX] (%i) -- Computing electrostatics serial\n", me);
+#endif
+
+        // compute energy+gradients in serial on rank 0 for full system
+
+        fix_mbx->mbxt_start(MBXT_ELE);
+        if (comm->me == 0) mbx_ele = ptr_mbx_full->Electrostatics(true);
+        fix_mbx->mbxt_stop(MBXT_ELE);
+        accumulate_f_full(true);
+    }
+
     if (mbx_parallel) {
 #ifdef _DEBUG
             printf("[MBX] (%i) -- Computing disp parallel\n", me);
@@ -229,29 +252,6 @@ void PairMBX::compute(int eflag, int vflag) {
 
     mbx_disp = mbx_disp_real + mbx_disp_pme;
 
-    if (mbx_parallel) {
-#ifdef _DEBUG
-        printf("[MBX] (%i) -- Computing electrostatics parallel\n", me);
-#endif
-
-        fix_mbx->mbxt_start(MBXT_ELE);
-        mbx_ele = ptr_mbx_local->ElectrostaticsMPIlocal(true, true);
-        fix_mbx->mbxt_stop(MBXT_ELE);
-        accumulate_f_local(true);
-
-    } else {
-#ifdef _DEBUG
-        printf("[MBX] (%i) -- Computing electrostatics serial\n", me);
-#endif
-
-        // compute energy+gradients in serial on rank 0 for full system
-
-        fix_mbx->mbxt_start(MBXT_ELE);
-        if (comm->me == 0) mbx_ele = ptr_mbx_full->Electrostatics(true);
-        fix_mbx->mbxt_stop(MBXT_ELE);
-        accumulate_f_full(true);
-    }
-
     mbx_total_energy = mbx_e1b + mbx_e2b + mbx_disp + mbx_buck + mbx_e3b + mbx_e4b + mbx_ele;
 
     for (int i = 0; i < 6; ++i) virial[i] += mbx_virial[i];

src/bblock/system.cpp

@@ -755,6 +755,53 @@ void System::SetPBC(std::vector<double> box) {
     SetPolfacs();
 }
 
+
+void System::SetPBCElectrostaticsMPIlocal(std::vector<double> box) {
+    // Check that the box has 0 or 9 components
+    if (box.size() != 9 && box.size() != 6 && box.size() != 0) {
+        std::string text = "Box size of " + std::to_string(box.size()) + " is not acceptable.";
+        throw CUException(__func__, __FILE__, __LINE__, text);
+    }
+
+#ifdef DEBUG
+    std::cerr << "Entered SetPBC():\n";
+    std::cerr << "Coordinate before fixing monomers:\n";
+    for (size_t i = 0; i < xyz_.size(); i++) {
+        std::cerr << xyz_[i] << " , ";
+    }
+    std::cerr << std::endl;
+#endif
+
+    // Set the box and the bool to use or not pbc
+    use_pbc_ = box.size();
+    if (use_pbc_) simcell_periodic_ = true;
+
+    box_ = box;
+    if (box.size() == 9) {
+        box_ = box;
+        box_ABCabc_ = BoxVecToBoxABCabc(box);
+    } else if (box.size() == 6) {
+        box_ABCabc_ = box;
+        box_ = BoxABCabcToBoxVec(box);
+    }
+
+    box_inverse_ = InvertUnitCell(box_);
+
+#ifdef DEBUG
+    std::cerr << "Coordinate after fixing monomers:\n";
+    for (size_t i = 0; i < xyz_.size(); i++) {
+        std::cerr << xyz_[i] << " , ";
+    }
+    std::cerr << std::endl;
+#endif
+
+    // Reset the virtual site positions, charges, pols and polfacs
+    SetVSites();
+    SetCharges();
+    SetPols();
+    SetPolfacs();
+}
+
 void System::SetXyz(std::vector<double> xyz) {
     // Make sure that the xyz of input has the right size
     if (xyz.size() != 3 * numsites_) {
@@ -3444,7 +3491,7 @@ double System::ElectrostaticsMPIlocal(bool do_grads, bool use_ghost) {
     if (islocal_.size() > 0) std::fill(grad_.begin(), grad_.end(), 0.0);
     std::fill(virial_.begin(), virial_.end(), 0.0);
 
-    SetPBC(box_);
+    SetPBCElectrostaticsMPIlocal(box_);
 
     energy_ = GetElectrostaticsMPIlocal(do_grads, use_ghost);
 

src/bblock/system.h

@@ -908,6 +908,19 @@ class System {
      */
     void SetPBC(std::vector<double> box = {});
 
+
+    /**
+     * Tells the system if we are in Periodic Boundary Conditions (PBC)
+     * or not. If the box is not passed as argument, it is set to
+     * a cubic box of 1000 Angstrom of side length
+     * @param[in] use_pbc Boolean setting if we use PBC (true) or not (false)
+     * @param[in] box Optional argument. Is a 9 component vector of double with
+     * the three main vectors of the cell: {v1x v1y v1z v2x v2y v2z v3x v3y v3z}
+     * 
+     * this function differs from SetPBC, because it does systools::FixMonomerCoordinates.
+     */
+    void SetPBCElectrostaticsMPIlocal(std::vector<double> box = {});
+
     /**
      * Sets the values for alpha, the grid density and the spline order when using PME to get the reciprocal space
      * contribution. Sets both the parameters for dispersion and electorstatics to be the same.