From 3d6ccab93b22b6541623956e763280a27b911882 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 17 Feb 2020 22:09:21 +0100 Subject: [PATCH 1/4] Fix periodic image bug The NameError was caused by a regression in 0a22eeb94. The `dim` variable was not defined, but the compiler warning was hidden by the wildcard import `from math import *`. --- src/python/espressomd/visualization_opengl.pyx | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/python/espressomd/visualization_opengl.pyx b/src/python/espressomd/visualization_opengl.pyx index be0b05b64cd..4b35a55d7bd 100644 --- a/src/python/espressomd/visualization_opengl.pyx +++ b/src/python/espressomd/visualization_opengl.pyx @@ -1220,9 +1220,9 @@ class openGLLive: for imz in range(-self.specs['periodic_images'][2], self.specs['periodic_images'][2] + 1): if imx != 0 or imy != 0 or imz != 0: - im = np.array([imx, imy, imz]) - draw_cylinder(self.particles['pos'][b[0]] + im * self.imPos[dim], self.particles['pos'][b[1]] + - im * self.imPos[dim], radius, col, mat, self.specs['quality_bonds']) + draw_cylinder(self.particles['pos'][b[0]] + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], + self.particles['pos'][b[1]] + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], + radius, col, mat, self.specs['quality_bonds']) else: # SPLIT BOND l = self.particles['pos'][b[0]] - self.particles['pos'][b[1]] @@ -1256,11 +1256,10 @@ class openGLLive: for imz in range(-self.specs['periodic_images'][2], self.specs['periodic_images'][2] + 1): if imx != 0 or imy != 0 or imz != 0: - im = np.array([imx, imy, imz]) - draw_cylinder(self.particles['pos'][b[0]] + im * self.imPos[dim], s0 + im * - self.imPos[dim], radius, col, mat, self.specs['quality_bonds']) - draw_cylinder(self.particles['pos'][b[1]] + im * self.imPos[dim], s1 + im * - self.imPos[dim], radius, col, mat, self.specs['quality_bonds']) + draw_cylinder(self.particles['pos'][b[0]] + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], + s0 + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], radius, col, mat, self.specs['quality_bonds']) + draw_cylinder(self.particles['pos'][b[1]] + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], + s1 + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], radius, col, mat, self.specs['quality_bonds']) def _redraw_sphere(self, pos, radius, quality): OpenGL.GL.glPushMatrix() From 4d555aa0ba8909e0bc2099c1a702cb8a9d87eec2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 17 Feb 2020 22:13:43 +0100 Subject: [PATCH 2/4] Remove wildcard import --- src/python/espressomd/visualization_opengl.pyx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/python/espressomd/visualization_opengl.pyx b/src/python/espressomd/visualization_opengl.pyx index 4b35a55d7bd..5dabe8ea1bc 100644 --- a/src/python/espressomd/visualization_opengl.pyx +++ b/src/python/espressomd/visualization_opengl.pyx @@ -17,7 +17,7 @@ import OpenGL.GLUT import OpenGL.GLU import OpenGL.GL -from math import * +import math import numpy as np import ctypes import os @@ -2020,7 +2020,7 @@ def rotation_helper(d): # the rotation axis is the cross product between z and d vz = np.cross([0.0, 0.0, 1.0], d) # get the angle using a dot product - angle = 180.0 / np.pi * acos(d[2] / np.linalg.norm(d)) + angle = 180.0 / np.pi * math.acos(d[2] / np.linalg.norm(d)) return angle, vz[0], vz[1] @@ -2217,7 +2217,7 @@ def draw_sphero_cylinder(posA, posB, radius, color, material, quality): if v == 0: ax = 57.2957795 else: - ax = 57.2957795 * acos(d[2] / v) + ax = 57.2957795 * math.acos(d[2] / v) if d[2] < 0.0: ax = -ax @@ -2574,8 +2574,8 @@ class Camera: return m def rotate_vector(self, vec, ang, axe): - sinhalf = sin(ang / 2) - coshalf = cos(ang / 2) + sinhalf = math.sin(ang / 2) + coshalf = math.cos(ang / 2) rx = axe[0] * sinhalf ry = axe[1] * sinhalf From fa36876add5a069a2d673ae7d57583e0a82a0ac7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 18 Feb 2020 17:12:22 +0100 Subject: [PATCH 3/4] Simplify code --- .../espressomd/visualization_opengl.pyx | 73 +++++++++---------- 1 file changed, 34 insertions(+), 39 deletions(-) diff --git a/src/python/espressomd/visualization_opengl.pyx b/src/python/espressomd/visualization_opengl.pyx index 5dabe8ea1bc..be7f474734d 100644 --- a/src/python/espressomd/visualization_opengl.pyx +++ b/src/python/espressomd/visualization_opengl.pyx @@ -413,6 +413,8 @@ class openGLLive: # HAS PERIODIC IMAGES self.has_images = any(i != 0 for i in self.specs['periodic_images']) + self.image_vectors = [ + list(range(-i, i + 1)) for i in self.specs['periodic_images']] # INITS self.system = system @@ -1141,15 +1143,13 @@ class openGLLive: self.particles['pos'][pid], radius, self.specs['quality_particles']) if self.has_images: - for imx in range(-self.specs['periodic_images'][0], - self.specs['periodic_images'][0] + 1): - for imy in range(-self.specs['periodic_images'][1], - self.specs['periodic_images'][1] + 1): - for imz in range(-self.specs['periodic_images'][2], - self.specs['periodic_images'][2] + 1): + for imx in self.image_vectors[0]: + for imy in self.image_vectors[1]: + for imz in self.image_vectors[2]: if imx != 0 or imy != 0 or imz != 0: - self._redraw_sphere( - self.particles['pos'][pid] + (imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2]), radius, self.specs['quality_particles']) + offset = [imx, imy, imz] * self.system.box_l + self._redraw_sphere(self.particles['pos'][pid] + offset, + radius, self.specs['quality_particles']) IF EXTERNAL_FORCES: if self.specs['ext_force_arrows'] or pid == self.dragId: @@ -1195,7 +1195,6 @@ class openGLLive: v, dtype=float) * sc, radius, col, self.materials['chrome'], self.specs['quality_arrows']) def _draw_bonds(self): - pIds = range(len(self.particles['pos'])) b2 = self.system.box_l[0] / 2.0 box_l2_sqr = pow(b2, 2.0) for b in self.bonds: @@ -1204,25 +1203,24 @@ class openGLLive: self.specs['bond_type_materials'], b[2])] radius = self._modulo_indexing( self.specs['bond_type_radius'], b[2]) - d = self.particles['pos'][b[0]] - self.particles['pos'][b[1]] - bondLen_sqr = d[0] * d[0] + d[1] * d[1] + d[2] * d[2] + xA = self.particles['pos'][b[0]] + xB = self.particles['pos'][b[1]] + dx = xB - xA + bondLen_sqr = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2] if bondLen_sqr < box_l2_sqr: # BOND COMPLETELY INSIDE BOX - draw_cylinder( - self.particles['pos'][ - b[0]], self.particles['pos'][b[1]], radius, - col, mat, self.specs['quality_bonds']) - for imx in range(-self.specs['periodic_images'][0], - self.specs['periodic_images'][0] + 1): - for imy in range(-self.specs['periodic_images'][1], - self.specs['periodic_images'][1] + 1): - for imz in range(-self.specs['periodic_images'][2], - self.specs['periodic_images'][2] + 1): - if imx != 0 or imy != 0 or imz != 0: - draw_cylinder(self.particles['pos'][b[0]] + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], - self.particles['pos'][b[1]] + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], - radius, col, mat, self.specs['quality_bonds']) + draw_cylinder(xA, xB, radius, col, mat, + self.specs['quality_bonds']) + if self.has_images: + for imx in self.image_vectors[0]: + for imy in self.image_vectors[1]: + for imz in self.image_vectors[2]: + if imx != 0 or imy != 0 or imz != 0: + offset = [imx, imy, imz] * \ + self.system.box_l + draw_cylinder(xA + offset, xB + offset, + radius, col, mat, self.specs['quality_bonds']) else: # SPLIT BOND l = self.particles['pos'][b[0]] - self.particles['pos'][b[1]] @@ -1249,17 +1247,17 @@ class openGLLive: draw_cylinder(self.particles['pos'][b[1]], s1, radius, col, mat, self.specs['quality_bonds']) - for imx in range(-self.specs['periodic_images'][0], - self.specs['periodic_images'][0] + 1): - for imy in range(-self.specs['periodic_images'][1], - self.specs['periodic_images'][1] + 1): - for imz in range(-self.specs['periodic_images'][2], - self.specs['periodic_images'][2] + 1): - if imx != 0 or imy != 0 or imz != 0: - draw_cylinder(self.particles['pos'][b[0]] + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], - s0 + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], radius, col, mat, self.specs['quality_bonds']) - draw_cylinder(self.particles['pos'][b[1]] + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], - s1 + imx * self.imPos[0] + imy * self.imPos[1] + imz * self.imPos[2], radius, col, mat, self.specs['quality_bonds']) + if self.has_images: + for imx in self.image_vectors[0]: + for imy in self.image_vectors[1]: + for imz in self.image_vectors[2]: + if imx != 0 or imy != 0 or imz != 0: + offset = [imx, imy, imz] * \ + self.system.box_l + draw_cylinder(xA + offset, s0 + offset, radius, col, mat, + self.specs['quality_bonds']) + draw_cylinder(xB + offset, s1 + offset, radius, col, mat, + self.specs['quality_bonds']) def _redraw_sphere(self, pos, radius, quality): OpenGL.GL.glPushMatrix() @@ -1618,9 +1616,6 @@ class openGLLive: self.depth = 0 - self.imPos = [np.array([self.system.box_l[0], 0, 0]), np.array( - [0, self.system.box_l[1], 0]), np.array([0, 0, self.system.box_l[2]])] - # LOOK FOR LB ACTOR if self.specs['LB_draw_velocity_plane'] or self.specs['LB_draw_nodes'] or self.specs['LB_draw_node_boundaries']: lb_types = [espressomd.lb.LBFluid] From ec59e22455c223173892c381d5bada4e72565c12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 18 Feb 2020 20:24:23 +0100 Subject: [PATCH 4/4] Fix bond connectivity issues at the box boundaries Calculate how to cut bonds that cross one side of the simulation box using an analytical solution instead of a heuristic algorithm. The heuristic was incorrectly configured, leading to wrong bond angles for cut bonds. This was particularly visible when drawing periodic images. Also, the old algorithm assumed the box was cubic, but now all directions are checked. --- .../espressomd/visualization_opengl.pyx | 108 ++++++++++-------- 1 file changed, 63 insertions(+), 45 deletions(-) diff --git a/src/python/espressomd/visualization_opengl.pyx b/src/python/espressomd/visualization_opengl.pyx index be7f474734d..a04c12b18df 100644 --- a/src/python/espressomd/visualization_opengl.pyx +++ b/src/python/espressomd/visualization_opengl.pyx @@ -1195,8 +1195,7 @@ class openGLLive: v, dtype=float) * sc, radius, col, self.materials['chrome'], self.specs['quality_arrows']) def _draw_bonds(self): - b2 = self.system.box_l[0] / 2.0 - box_l2_sqr = pow(b2, 2.0) + box_l_2 = self.system.box_l / 2.0 for b in self.bonds: col = self._modulo_indexing(self.specs['bond_type_colors'], b[2]) mat = self.materials[self._modulo_indexing( @@ -1206,9 +1205,9 @@ class openGLLive: xA = self.particles['pos'][b[0]] xB = self.particles['pos'][b[1]] dx = xB - xA - bondLen_sqr = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2] - if bondLen_sqr < box_l2_sqr: + if abs(dx[0]) < box_l_2[0] and abs( + dx[1]) < box_l_2[1] and abs(dx[2]) < box_l_2[2]: # BOND COMPLETELY INSIDE BOX draw_cylinder(xA, xB, radius, col, mat, self.specs['quality_bonds']) @@ -1222,42 +1221,29 @@ class openGLLive: draw_cylinder(xA + offset, xB + offset, radius, col, mat, self.specs['quality_bonds']) else: - # SPLIT BOND - l = self.particles['pos'][b[0]] - self.particles['pos'][b[1]] - l0 = self.particles['pos'][b[0]] - hits = 0 - for i in range(6): - lineBoxNDot = float(np.dot(l, self.box_n[i])) - if lineBoxNDot == 0: - continue - s = l0 - \ - np.dot(l0 - self.box_p[i], - self.box_n[i]) / lineBoxNDot * l - if self._is_inside_box(s): - if lineBoxNDot < 0: - s0 = s - else: - s1 = s - hits += 1 - if hits >= 2: - break - if hits >= 2: - draw_cylinder(self.particles['pos'][b[0]], s0, radius, col, - mat, self.specs['quality_bonds']) - draw_cylinder(self.particles['pos'][b[1]], s1, radius, col, - mat, self.specs['quality_bonds']) - - if self.has_images: - for imx in self.image_vectors[0]: - for imy in self.image_vectors[1]: - for imz in self.image_vectors[2]: - if imx != 0 or imy != 0 or imz != 0: - offset = [imx, imy, imz] * \ - self.system.box_l - draw_cylinder(xA + offset, s0 + offset, radius, col, mat, - self.specs['quality_bonds']) - draw_cylinder(xB + offset, s1 + offset, radius, col, mat, - self.specs['quality_bonds']) + # BOND CROSSES THE BOX BOUNDARIES + d = self._cut_bond(xA, dx) + if d is np.inf: + continue + + sA = xA + d * dx + sB = xB - (1 - d) * dx + draw_cylinder(xA, sA, radius, col, mat, + self.specs['quality_bonds']) + draw_cylinder(xB, sB, radius, col, mat, + self.specs['quality_bonds']) + + if self.has_images: + for imx in self.image_vectors[0]: + for imy in self.image_vectors[1]: + for imz in self.image_vectors[2]: + if imx != 0 or imy != 0 or imz != 0: + offset = [imx, imy, imz] * \ + self.system.box_l + draw_cylinder(xA + offset, sA + offset, radius, col, mat, + self.specs['quality_bonds']) + draw_cylinder(xB + offset, sB + offset, radius, col, mat, + self.specs['quality_bonds']) def _redraw_sphere(self, pos, radius, quality): OpenGL.GL.glPushMatrix() @@ -1266,12 +1252,44 @@ class openGLLive: OpenGL.GL.OpenGL.GL.glPopMatrix() # HELPER TO DRAW PERIODIC BONDS - def _is_inside_box(self, p): - eps = 1e-5 + def _cut_bond(self, xA, dx): + """ + Algorithm for the line-plane intersection problem. Given the unfolded + positions of two particles, determine: 1) the box image that minimizes + the folded bond length, 2) which side of the simulation box is crossed + by the bond, 3) how much of the bond is inside the unit cell starting + from the first particle. That scalar is returned by the function, and + can be used to calculate the coordinates of the line-plane + intersection point. The algorithm can be found at + https://en.wikipedia.org/w/index.php?title=Line%E2%80%93plane_intersection&oldid=940427752#Algebraic_form + """ + # corner case: the two particles occupy the same position + if np.dot(dx, dx) < 1e-9: + return np.inf + # find the box image that minimizes the unfolded bond length + shift = np.rint(dx / self.system.box_l) + # get the unfolded bond vector + dx -= shift * self.system.box_l + # find which side of the simulation box is crossed by the bond + best_d = np.inf + best_i = 0 for i in range(3): - if p[i] < -eps or p[i] > eps + self.last_box_l[i]: - return False - return True + if dx[i] == 0: + continue # corner case: bond is parallel to a face + elif dx[i] > 0: + p0_i = self.system.box_l[i] + else: + p0_i = 0 + # calculate the length of the bond that is inside the box using + # an optimized version of `np.dot(p0 - x0, n) / np.dot(dx, n)` + # where the dot products decay to an array access, because `n` + # is always a unit vector in rectangular boxes and its sign + # is canceled out by the division + d = (p0_i - xA[i]) / dx[i] + if d < best_d: + best_d = d + best_i = i + return best_d # ARROWS IN A PLANE FOR LB VELOCITIES def _draw_lb_vel(self):