Merge #3128

3128: Icc vs energy fix r=RudolfWeeber a=itischler Fixes #3079 Description of changes: - changes in virtual sites and icc are now applied in on_observable_calc PR Checklist ------------ - [ ] Tests? - [ ] Interface - [ ] Core - [ ] Docs? Co-authored-by: itischler <[email protected]>
espressomd · Sep 9, 2019 · b814c25 · b814c25
2 parents 4d69cd2 + 207ba31
commit b814c25
Show file tree

Hide file tree

Showing 5 changed files with 33 additions and 1 deletion.
diff --git a/src/core/communication.cpp b/src/core/communication.cpp
@@ -612,6 +612,7 @@ void mpi_iccp3m_init_slave(const iccp3m_struct &iccp3m_cfg_) {
 #ifdef ELECTROSTATICS
   iccp3m_cfg = iccp3m_cfg_;
 
+  on_particle_charge_change();
   check_runtime_errors(comm_cart);
 #endif
 }
@@ -622,6 +623,7 @@ int mpi_iccp3m_init() {
 #ifdef ELECTROSTATICS
   mpi_call(mpi_iccp3m_init_slave, iccp3m_cfg);
 
+  on_particle_charge_change();
   return check_runtime_errors(comm_cart);
 #else
   return 0;

diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp
@@ -9,6 +9,7 @@ double coulomb_cutoff;
 #include "debug.hpp"
 #include "electrostatics_magnetostatics/debye_hueckel.hpp"
 #include "electrostatics_magnetostatics/elc.hpp"
+#include "electrostatics_magnetostatics/icc.hpp"
 #include "electrostatics_magnetostatics/mmm1d.hpp"
 #include "electrostatics_magnetostatics/mmm2d.hpp"
 #include "electrostatics_magnetostatics/p3m.hpp"
@@ -181,6 +182,8 @@ void integrate_sanity_check() {
 }
 
 void on_observable_calc() {
+  iccp3m_iteration(local_cells.particles(),
+                   cell_structure.ghost_cells().particles());
   switch (coulomb.method) {
 #ifdef P3M
   case COULOMB_ELC_P3M:

diff --git a/src/core/event.cpp b/src/core/event.cpp
@@ -178,7 +178,15 @@ void on_observable_calc() {
    * information */
 
   cells_update_ghosts();
-
+  if (recalc_forces) {
+#ifdef VIRTUAL_SITES
+    if (virtual_sites()->is_relative()) {
+      ghost_communicator(&cell_structure.update_ghost_pos_comm);
+    }
+    virtual_sites()->update();
+#endif
+    cells_update_ghosts();
+  }
 #ifdef ELECTROSTATICS
   if (reinit_electrostatics) {
     EVENT_TRACE(fprintf(stderr, "%d: reinit_electrostatics\n", this_node));

diff --git a/testsuite/python/icc.py b/testsuite/python/icc.py
@@ -101,6 +101,16 @@ def runTest(self):
         # Result
         self.assertAlmostEqual(1, induced_dipole / testcharge_dipole, places=4)
 
+        # Test applying changes
+        enegry_pre_change = S.analysis.energy()['total']
+        pressure_pre_change = S.analysis.pressure()['total']
+        icc.set_params(sigmas=[2.0] * nicc_tot)
+        icc.set_params(epsilons=[20.0] * nicc_tot)
+        enegry_post_change = S.analysis.energy()['total']
+        pressure_post_change = S.analysis.pressure()['total']
+        self.assertNotAlmostEqual(enegry_pre_change, enegry_post_change)
+        self.assertNotAlmostEqual(pressure_pre_change, pressure_post_change)
+
 
 if __name__ == "__main__":
     ut.main()
diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py
@@ -290,6 +290,15 @@ def run_test_lj(self):
             verify_lj_forces(system, 1E-10, 3 *
                              np.arange(int(n / 2), dtype=int))
 
+        # Test applying changes
+        enegry_pre_change = system.analysis.energy()['total']
+        pressure_pre_change = system.analysis.pressure()['total']
+        system.part[0].pos = system.part[0].pos + (2.2, -1.4, 4.2)
+        enegry_post_change = system.analysis.energy()['total']
+        pressure_post_change = system.analysis.pressure()['total']
+        self.assertNotAlmostEqual(enegry_pre_change, enegry_post_change)
+        self.assertNotAlmostEqual(pressure_pre_change, pressure_post_change)
+
         # Turn off lj interaction
         system.non_bonded_inter[0, 0].lennard_jones.set_params(
             epsilon=0, sigma=0, cutoff=0, shift=0)