Combination of fixing rotation, dihedral and virtual sites broken with steepest descent algorithm

[KaiSzuttor]

The following script shows the bug:

import sys

import numpy as np

import espressomd

system = espressomd.System(box_l=[10.0] * 3)
system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
system.time_step = 0.01
system.cell_system.skin = 0.4

positions = np.array([[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0], [1.0, 2.0, 1.0]])
system.part.add(pos=positions)

dihedral = espressomd.interactions.Dihedral(bend=10.0, mult=1, phase=np.pi/4.0)
system.bonded_inter.add(dihedral)
system.part[1].add_bond((dihedral, 0, 2, 3))
system.part[1].rotation = [False] * 3 # <------ this line breaks the dihedral bond

system.part[0].vs_auto_relate_to(1)
system.part[2].vs_auto_relate_to(1)
system.part[3].vs_auto_relate_to(1)

system.integrator.set_steepest_descent(f_max=0, gamma=20,
                                       max_displacement=0.01)
energy = system.analysis.energy()['total']
system.integrator.run(1)

[KaiSzuttor]

if the marked line is changed to system.part[1].rotation = [True] * 3, the bond does not break.

[RudolfWeeber]

Please try the floolowing:

In rotation.hpp:local_rotate_particle or local_rotate_particle_body, the rotation axis in body fied frame needs to be masked by the rotation flags fo the particle.

Use mask(p.p.rotation,axis_in_body_frame)

[KaiSzuttor]

it's already masked in line 111 of rotation.hpp

[KaiSzuttor]

virtual sites?

Am 10.11.2020 um 20:45 schrieb Florian Weik notifications@github.com:

 In your example you have nothing keeping the particles of the dihedral together. Why shouldn't it break?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[jngrad]

With rotation disabled, the real particles are collapsed onto the real particle.

diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp
index 9bac61bfc..aa5ea496d 100644
--- a/src/core/integrate.cpp
+++ b/src/core/integrate.cpp
@@ -200,2 +200,3 @@ int integrate(int n_steps, int reuse_forces) {
     ESPRESSO_PROFILER_MARK_BEGIN("Initial Force Calculation");
+    std::cout << "recalc_forces\n";
     lb_lbcoupling_deactivate();
@@ -235,2 +236,3 @@ int integrate(int n_steps, int reuse_forces) {
   for (int step = 0; step < n_steps; step++) {
+    std::cout << "integrate\n";
     ESPRESSO_PROFILER_CXX_MARK_LOOP_ITERATION(integration_loop, step);
diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp
index 293f1ec84..fe5470bbb 100644
--- a/src/core/virtual_sites/VirtualSitesRelative.cpp
+++ b/src/core/virtual_sites/VirtualSitesRelative.cpp
@@ -151,4 +151,9 @@ auto constraint_stress(
 void VirtualSitesRelative::update() const {
+  std::cout << "VirtualSitesRelative::update()\n";
   cell_structure.ghosts_update(Cells::DATA_PART_POSITION |
                                Cells::DATA_PART_MOMENTUM);
+  for (auto &p : cell_structure.local_particles()) {
+    printf("original %d : %.1f %.1f %.1f", p.p.identity, p.r.p[0], p.r.p[1], p.r.p[2]);
+    if (p.p.is_virtual) printf(" (virtual)\n"); else printf(" (real)\n");
+  }

@@ -174,2 +179,6 @@ void VirtualSitesRelative::update() const {
   } // namespace
+  for (auto &p : cell_structure.local_particles()) {
+    printf("shift    %d : %.1f %.1f %.1f", p.p.identity, p.r.p[0], p.r.p[1], p.r.p[2]);
+    if (p.p.is_virtual) printf(" (virtual)\n"); else printf(" (real)\n");
+  }
 }

output:

[...]
integrate
VirtualSitesRelative::update()
original 0 : 0.0 0.0 0.0 (virtual)
original 1 : 0.0 1.0 0.0 (real)
original 2 : 1.0 1.0 0.0 (virtual)
original 3 : 1.0 2.0 1.0 (virtual)
shift    0 : 0.0 1.0 0.0 (virtual)
shift    1 : 0.0 1.0 0.0 (real)
shift    2 : 0.0 1.0 0.0 (virtual)
shift    3 : 0.0 1.0 0.0 (virtual)
VirtualSitesRelative::update()
original 0 : 0.0 1.0 0.0 (virtual)
original 1 : 0.0 1.0 0.0 (real)
original 2 : 0.0 1.0 0.0 (virtual)
original 3 : 0.0 1.0 0.0 (virtual)
shift    0 : 0.0 1.0 0.0 (virtual)
shift    1 : 0.0 1.0 0.0 (real)
shift    2 : 0.0 1.0 0.0 (virtual)
shift    3 : 0.0 1.0 0.0 (virtual)

[jngrad]

This bug got fixed when the hand-written quaternion code was replaced by the boost quaternion code in f0f72e32bf7d2e0d38b6d31a62edd7840446620c in 4.2.0.