Merge #3216 #3239 #3253 #3258

3216: Refactor ghosts.cpp r=KaiSzuttor a=hirschsn More refactoring that builds upon #3212 Description of *major* changes: - remove GhostCommunication::mpi_comm as it is not used in ghosts.cpp, - make GhostCommunication::part_lists a std::vector, - remove static variables s_buffer and r_buffer, - factor out memory handling, - change loops to range based for, - use boost::mpi. - Replace the manual poststore and prefetch loops by find_ifs Mainly, ghosts.cpp now defines CommBuf, which is a container for the data to be sent or received as well as two classes (Archiver and BondArchiver), that insert and extract the memory from CommBuf. Some of these changes, like the removal of static variables, is necessary for my implementation of asynchronous ghost communication. PR Checklist ------------ - [ ] Tests? - [ ] Interface - [ ] Core - [ ] Docs? 3239: Added test criteria for the charged_system-2 tutorial r=RudolfWeeber,jngrad a=reinaual 3253: Refactor NpT public interface r=fweik a=jngrad Description of changes: - remove the silent conversion of the incorrect input parameter `dimension=[0,0,0]` to `[1,1,1]` in the core (bypassing sanity checks), now the checks will throw an exception for fixed-volume NpT; the original behavior was counter-intuitive and undocumented until 2 days ago - remove the automatic decay of NpT to NVT upon initialization of NpT with incorrect parameters - remove unused `p_inst_av` variable (average instantaneous pressure) - cleanup integrator documentation 3258: CMake minor fixes r=fweik a=jngrad Description of changes: - change next milestone to 4.2 - load `GNUInstallDirs` to make standard GNU paths accessible from CMake variables - simplify CMake logic and install in `python3.X` folder instead of the deprecated `python3` folder - add extra check to make sure install paths are correctly configured (all python and shared object files must be inside the package `espressomd`) Co-authored-by: Steffen Hirschmann <[email protected]> Co-authored-by: Florian Weik <[email protected]> Co-authored-by: Alexander Reinauer <[email protected]> Co-authored-by: Jean-Noël Grad <[email protected]>
espressomd · Oct 18, 2019 · 628b15a · 628b15a
5 parents 192b788 + f21d717 + 129029a + 98e2749 + db8cb80
commit 628b15a
Show file tree

Hide file tree

Showing 31 changed files with 583 additions and 617 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -20,6 +20,7 @@
 
 cmake_minimum_required(VERSION 3.4)
 include(FeatureSummary)
+include(GNUInstallDirs)
 project(ESPResSo)
 include(cmake/FindPythonModule.cmake)
 if(NOT ${CMAKE_VERSION} VERSION_LESS "3.12")
@@ -28,7 +29,7 @@ endif()
 
 enable_language(CXX)
 
-set(PROJECT_VERSION "5.0-dev")
+set(PROJECT_VERSION "4.2-dev")
 
 #
 # CMake internal vars
@@ -241,7 +242,7 @@ if(WITH_PYTHON)
     execute_process(
       COMMAND
         ${PYTHON_EXECUTABLE} -c
-        "import distutils.sysconfig as cg; print(cg.get_python_lib(1,0,prefix='${CMAKE_INSTALL_EXEC_PREFIX}'))"
+        "import distutils.sysconfig as cg; print(cg.get_python_lib(prefix='${CMAKE_INSTALL_PREFIX}', plat_specific=True, standard_lib=False).replace('${CMAKE_INSTALL_PREFIX}/', '', 1))"
       OUTPUT_VARIABLE PYTHON_INSTDIR
       OUTPUT_STRIP_TRAILING_WHITESPACE)
   endif(NOT PYTHON_INSTDIR)
@@ -398,10 +399,6 @@ endif()
 # Paths
 #
 
-if(NOT DEFINED BINDIR)
-  set(BINDIR "bin")
-endif(NOT DEFINED BINDIR)
-
 set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/${PYTHON_INSTDIR}/espressomd")
 
 #

diff --git a/doc/sphinx/analysis.rst b/doc/sphinx/analysis.rst
@@ -669,7 +669,7 @@ which the observable should take into account.
 
 The current value of an observable can be obtained using its calculate()-method::
 
-    print(par_pos.calculate())
+    print(part_pos.calculate())
 
 .. _Available observables:
 

diff --git a/doc/tutorials/02-charged_system/02-charged_system-1.ipynb b/doc/tutorials/02-charged_system/02-charged_system-1.ipynb
@@ -19,7 +19,7 @@
     "#define EXTERNAL_FORCES\n",
     "#define MASS\n",
     "#define ELECTROSTATICS\n",
-    "#define LENNARD_JONES\n",
+    "#define WCA\n",
     "```"
    ]
   },
@@ -46,7 +46,7 @@
     "plt.ion()\n",
     "\n",
     "# Print enabled features\n",
-    "required_features = [\"EXTERNAL_FORCES\", \"MASS\", \"ELECTROSTATICS\", \"LENNARD_JONES\"]\n",
+    "required_features = [\"EXTERNAL_FORCES\", \"MASS\", \"ELECTROSTATICS\", \"WCA\"]\n",
     "espressomd.assert_features(required_features)\n",
     "print(espressomd.features())\n",
     "\n",
@@ -67,12 +67,8 @@
     "types       = {\"Anion\":          0, \"Cation\": 1}\n",
     "numbers     = {\"Anion\": n_ionpairs, \"Cation\": n_ionpairs}\n",
     "charges     = {\"Anion\":       -1.0, \"Cation\": 1.0}\n",
-    "lj_sigmas   = {\"Anion\":        1.0, \"Cation\": 1.0}\n",
-    "lj_epsilons = {\"Anion\":        1.0, \"Cation\": 1.0}\n",
-    "\n",
-    "WCA_cut = 2.**(1. / 6.)\n",
-    "lj_cuts     = {\"Anion\":  WCA_cut * lj_sigmas[\"Anion\"],\n",
-    "               \"Cation\": WCA_cut * lj_sigmas[\"Cation\"]}"
+    "wca_sigmas   = {\"Anion\":        1.0, \"Cation\": 1.0}\n",
+    "wca_epsilons = {\"Anion\":        1.0, \"Cation\": 1.0}"
    ]
   },
   {
@@ -139,12 +135,11 @@
    "metadata": {},
    "source": [
     "Before we can really start the simulation, we have to specify the\n",
-    "interactions between our particles. We already defined the Lennard-Jones parameters at the beginning,\n",
+    "interactions between our particles. We already defined the WCA parameters at the beginning,\n",
     "what is left is to specify the combination rule and to iterate over particle type pairs. For simplicity, \n",
     "we implement only the *Lorentz-Berthelot* rules. \n",
     "We pass our interaction pair to <tt>system.non_bonded_inter[\\*,\\*]</tt> and set the \n",
-    "pre-calculated LJ parameters <tt>epsilon</tt>, <tt>sigma</tt> and <tt>cutoff</tt>. With <tt>shift=\"auto\"</tt>,\n",
-    "we shift the interaction potential to the cutoff so that $U_\\mathrm{LJ}(r_\\mathrm{cutoff})=0$."
+    "pre-calculated WCA parameters <tt>epsilon</tt> and <tt>sigma</tt>."
    ]
   },
   {
@@ -165,14 +160,13 @@
     "    else:\n",
     "        return ValueError(\"No combination rule defined\")\n",
     "\n",
-    "# Lennard-Jones interactions parameters\n",
+    "# WCA interactions parameters\n",
     "for s in [[\"Anion\", \"Cation\"], [\"Anion\", \"Anion\"], [\"Cation\", \"Cation\"]]:\n",
-    "    lj_sig = combination_rule_sigma(\"Berthelot\", lj_sigmas[s[0]], lj_sigmas[s[1]])\n",
-    "    lj_cut = combination_rule_sigma(\"Berthelot\", lj_cuts[s[0]], lj_cuts[s[1]])\n",
-    "    lj_eps = combination_rule_epsilon(\"Lorentz\", lj_epsilons[s[0]], lj_epsilons[s[1]])\n",
+    "    wca_sig = combination_rule_sigma(\"Berthelot\", wca_sigmas[s[0]], wca_sigmas[s[1]])\n",
+    "    wca_eps = combination_rule_epsilon(\"Lorentz\", wca_epsilons[s[0]], wca_epsilons[s[1]])\n",
     "\n",
-    "    system.non_bonded_inter[types[s[0]], types[s[1]]].lennard_jones.set_params(\n",
-    "        epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift=\"auto\")"
+    "    system.non_bonded_inter[types[s[0]], types[s[1]]].wca.set_params(\n",
+    "        epsilon=wca_eps, sigma=wca_sig)"
    ]
   },
   {
@@ -182,7 +176,7 @@
     "## 3 Equilibration\n",
     "\n",
     "With randomly positioned particles, we most likely have huge overlap and the strong repulsion will\n",
-    "cause the simulation to crash. The next step in our script therefore is a suitable LJ equilibration.\n",
+    "cause the simulation to crash. The next step in our script therefore is a suitable WCA equilibration.\n",
     "This is known to be a tricky part of a simulation and several approaches exist to reduce the particle overlap.\n",
     "Here, we use the steepest descent algorithm and cap the maximal particle displacement per integration step\n",
     "to 1% of sigma.\n",
@@ -197,8 +191,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Lennard-Jones Equilibration\n",
-    "max_sigma = max(lj_sigmas.values())\n",
+    "# WCA Equilibration\n",
+    "max_sigma = max(wca_sigmas.values())\n",
     "min_dist = 0.0\n",
     "system.minimize_energy.init(f_max=0, gamma=10.0, max_steps=10,\n",
     "                            max_displacement=max_sigma * 0.01)\n",
@@ -423,7 +417,7 @@
     "\n",
     "To make your life easier, don't try to equilibrate randomly positioned particles,\n",
     "but set them up in a crystal structure close to equilibrium. If you do it right,\n",
-    "you don't even need the Lennard-Jones equilibration. \n",
+    "you don't even need the WCA equilibration. \n",
     "To speed things up, don't go further than 1000 particles and use a P$^3$M accuracy of $10^{-2}$.\n",
     "Your RDF should look like the plot in figure 2. If you get stuck,\n",
     "you can look at the solution script <tt>/doc/tutorials/02-charged_system/scripts/nacl_units.py</tt> (or <tt>nacl_units_vis.py</tt> with visualization)."

diff --git a/doc/tutorials/02-charged_system/scripts/nacl.py b/doc/tutorials/02-charged_system/scripts/nacl.py
@@ -22,7 +22,7 @@
 from espressomd import assert_features, electrostatics
 import numpy
 
-assert_features(["ELECTROSTATICS", "LENNARD_JONES"])
+assert_features(["ELECTROSTATICS", "WCA"])
 
 print("\n--->Setup system")
 
@@ -43,12 +43,8 @@
 types = {"Anion": 0, "Cation": 1}
 numbers = {"Anion": n_ionpairs, "Cation": n_ionpairs}
 charges = {"Anion": -1.0, "Cation": 1.0}
-lj_sigmas = {"Anion": 1.0, "Cation": 1.0}
-lj_epsilons = {"Anion": 1.0, "Cation": 1.0}
-
-WCA_cut = 2.**(1. / 6.)
-lj_cuts = {"Anion": WCA_cut * lj_sigmas["Anion"],
-           "Cation": WCA_cut * lj_sigmas["Cation"]}
+wca_sigmas = {"Anion": 1.0, "Cation": 1.0}
+wca_epsilons = {"Anion": 1.0, "Cation": 1.0}
 
 # Setup System
 box_l = (n_part / density)**(1. / 3.)
@@ -83,18 +79,17 @@ def combination_rule_sigma(rule, sig1, sig2):
 
 # Lennard-Jones interactions parameters
 for s in [["Anion", "Cation"], ["Anion", "Anion"], ["Cation", "Cation"]]:
-    lj_sig = combination_rule_sigma(
-        "Berthelot", lj_sigmas[s[0]], lj_sigmas[s[1]])
-    lj_cut = combination_rule_sigma("Berthelot", lj_cuts[s[0]], lj_cuts[s[1]])
-    lj_eps = combination_rule_epsilon(
-        "Lorentz", lj_epsilons[s[0]], lj_epsilons[s[1]])
+    wca_sig = combination_rule_sigma(
+        "Berthelot", wca_sigmas[s[0]], wca_sigmas[s[1]])
+    wca_eps = combination_rule_epsilon(
+        "Lorentz", wca_epsilons[s[0]], wca_epsilons[s[1]])
 
-    system.non_bonded_inter[types[s[0]], types[s[1]]].lennard_jones.set_params(
-        epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
+    system.non_bonded_inter[types[s[0]], types[s[1]]].wca.set_params(
+        epsilon=wca_eps, sigma=wca_sig)
 
 
-print("\n--->Lennard-Jones Equilibration")
-max_sigma = max(lj_sigmas.values())
+print("\n--->WCA Equilibration")
+max_sigma = max(wca_sigmas.values())
 min_dist = 0.0
 system.minimize_energy.init(f_max=0, gamma=10, max_steps=10,
                             max_displacement=max_sigma * 0.01)

diff --git a/src/core/communication.cpp b/src/core/communication.cpp
@@ -537,7 +537,7 @@ void mpi_bcast_nptiso_geom() {
 void mpi_bcast_nptiso_geom_slave(int, int) {
   MPI_Bcast(&nptiso.geometry, 1, MPI_INT, 0, comm_cart);
   MPI_Bcast(&nptiso.dimension, 1, MPI_INT, 0, comm_cart);
-  MPI_Bcast(&nptiso.cubic_box, 1, MPI_INT, 0, comm_cart);
+  MPI_Bcast(&nptiso.cubic_box, 1, MPI_LOGICAL, 0, comm_cart);
   MPI_Bcast(&nptiso.non_const_dim, 1, MPI_INT, 0, comm_cart);
 }
 

diff --git a/src/core/domain_decomposition.cpp b/src/core/domain_decomposition.cpp
@@ -281,9 +281,7 @@ void dd_prepare_comm(GhostCommunicator *comm, int data_parts,
           comm->comm[cnt].node = this_node;
 
           /* Buffer has to contain Send and Recv cells -> factor 2 */
-          comm->comm[cnt].part_lists =
-              (Cell **)Utils::malloc(2 * n_comm_cells[dir] * sizeof(Cell *));
-          comm->comm[cnt].n_part_lists = 2 * n_comm_cells[dir];
+          comm->comm[cnt].part_lists.resize(2 * n_comm_cells[dir]);
           /* prepare folding of ghost positions */
           if ((data_parts & GHOSTTRANS_POSSHFTD) &&
               local_geo.boundary()[2 * dir + lr] != 0) {
@@ -294,7 +292,7 @@ void dd_prepare_comm(GhostCommunicator *comm, int data_parts,
           /* fill send comm cells */
           lc[dir] = hc[dir] = 1 + lr * (dd.cell_grid[dir] - 1);
 
-          dd_fill_comm_cell_lists(comm->comm[cnt].part_lists, lc, hc);
+          dd_fill_comm_cell_lists(comm->comm[cnt].part_lists.data(), lc, hc);
 
           /* fill recv comm cells */
           lc[dir] = hc[dir] = 0 + (1 - lr) * (dd.cell_grid[dir] + 1);
@@ -313,9 +311,7 @@ void dd_prepare_comm(GhostCommunicator *comm, int data_parts,
             if ((node_pos[dir] + i) % 2 == 0) {
               comm->comm[cnt].type = GHOST_SEND;
               comm->comm[cnt].node = node_neighbors[2 * dir + lr];
-              comm->comm[cnt].part_lists =
-                  (Cell **)Utils::malloc(n_comm_cells[dir] * sizeof(Cell *));
-              comm->comm[cnt].n_part_lists = n_comm_cells[dir];
+              comm->comm[cnt].part_lists.resize(n_comm_cells[dir]);
               /* prepare folding of ghost positions */
               if ((data_parts & GHOSTTRANS_POSSHFTD) &&
                   local_geo.boundary()[2 * dir + lr] != 0) {
@@ -325,21 +321,21 @@ void dd_prepare_comm(GhostCommunicator *comm, int data_parts,
 
               lc[dir] = hc[dir] = 1 + lr * (dd.cell_grid[dir] - 1);
 
-              dd_fill_comm_cell_lists(comm->comm[cnt].part_lists, lc, hc);
+              dd_fill_comm_cell_lists(comm->comm[cnt].part_lists.data(), lc,
+                                      hc);
               cnt++;
             }
           if (box_geo.periodic(dir) ||
               (local_geo.boundary()[2 * dir + (1 - lr)] == 0))
             if ((node_pos[dir] + (1 - i)) % 2 == 0) {
               comm->comm[cnt].type = GHOST_RECV;
               comm->comm[cnt].node = node_neighbors[2 * dir + (1 - lr)];
-              comm->comm[cnt].part_lists =
-                  (Cell **)Utils::malloc(n_comm_cells[dir] * sizeof(Cell *));
-              comm->comm[cnt].n_part_lists = n_comm_cells[dir];
+              comm->comm[cnt].part_lists.resize(n_comm_cells[dir]);
 
               lc[dir] = hc[dir] = (1 - lr) * (dd.cell_grid[dir] + 1);
 
-              dd_fill_comm_cell_lists(comm->comm[cnt].part_lists, lc, hc);
+              dd_fill_comm_cell_lists(comm->comm[cnt].part_lists.data(), lc,
+                                      hc);
               cnt++;
             }
         }
@@ -370,7 +366,7 @@ void dd_revert_comm_order(GhostCommunicator *comm) {
     else if (comm->comm[i].type == GHOST_RECV)
       comm->comm[i].type = GHOST_SEND;
     else if (comm->comm[i].type == GHOST_LOCL) {
-      nlist2 = comm->comm[i].n_part_lists / 2;
+      nlist2 = comm->comm[i].part_lists.size() / 2;
       for (j = 0; j < nlist2; j++) {
         auto tmplist = comm->comm[i].part_lists[j];
         comm->comm[i].part_lists[j] = comm->comm[i].part_lists[j + nlist2];

diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp
@@ -334,7 +334,7 @@ void calc_long_range_force(const ParticleRange &particles) {
     if (this_node == 0) {
       p3m_gpu_add_farfield_force();
     }
-    /* there is no NPT handling here as long as we cannot compute energies.g
+    /* there is no NPT handling here as long as we cannot compute energies.
        This is checked in integrator_npt_sanity_checks() when integration
        starts. */
     break;