From 4ef3e307559d8404396b837ccf3ec61caa1f6931 Mon Sep 17 00:00:00 2001 From: "dr.yuri" Date: Mon, 27 May 2019 10:54:44 +0300 Subject: [PATCH] lattice translation with S.Harper basis --- src/lattice.cpp | 251 ++++++++++++++++++++++++------------------------ 1 file changed, 127 insertions(+), 124 deletions(-) diff --git a/src/lattice.cpp b/src/lattice.cpp index af025033ac..5a3fea58a2 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -694,134 +694,137 @@ std::pair> HexLattice::distance(Position r, Direction u, const std::array& i_xyz) const { - double cosa,sina,u_xy; - // Index + 1 of translation direction - int i_trans {1}; - xt :: xarray translation_matrix; - if (hextype==0){ - cosa = std::sqrt(3.0) / 2.0; - sina = 0.5; - u_xy = u.y; - //Translation matrix 4x3 {{beta},{gamma},{oy},{oz}} for OY case - translation_matrix = {{1, 0, 0}, - {1, -1, 0}, - {0, 1, 0}, - {0, 0, 1}}; - } else { - cosa = 0.5; - sina = std::sqrt(3.0) / 2.0; - u_xy = u.x; - //Translation matrix 4x3 {{beta},{gamma},{ox},{oz}} for OX case - translation_matrix = {{0, 1, 0}, - {1, -1, 0}, - {1, 0, 0}, - {0, 0, 1}}; - } - double beta_dir = u.x * cosa + u.y * sina; - double gamma_dir = u.x * cosa - u.y * sina; + //Short description + //OY - orientation: + // basis0 = (1, 0) + // basis1 = (-1/sqrt(3), 1) = +120 degrees from basis0 + // beta = (sqrt(3)/2, 1/2) = +30 degrees from basis0 + // gamma = (sqrt(3)/2, -1/2) = -60 degrees from beta + // delta = (0, 1) = +60 degrees from beta + //OX - orientation: + // basis0 = (1/sqrt(3), -1) + // basis1 = (0, 1) = +120 degrees from basis0 + // beta = (1, 0) = +30 degrees from basis0 + // gamma = (1/2, -sqrt(3)/2) = -60 degrees from beta + // delta = (1/2, sqrt(3)/2) = +60 degrees from beta + //OZ be considered separetly + double beta_dir; + double gamma_dir; + double delta_dir; + if (hextype==0){ + beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; + gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0; + delta_dir = u.y; + } + else{ + beta_dir = u.x; + gamma_dir = u.x / 2.0 - u.y * std::sqrt(3.0) / 2.0; + delta_dir = u.x / 2.0 + u.y * std::sqrt(3.0) / 2.0; + } - // Note that hexagonal lattice distance calculations are performed - // using the particle's coordinates relative to the neighbor lattice - // cells, not relative to the particle's current cell. This is done - // because there is significant disagreement between neighboring cells - // on where the lattice boundary is due to finite precision issues. + // Note that hexagonal lattice distance calculations are performed + // using the particle's coordinates relative to the neighbor lattice + // cells, not relative to the particle's current cell. This is done + // because there is significant disagreement between neighboring cells + // on where the lattice boundary is due to finite precision issues. - //beta direction - double d {INFTY}; - std::array lattice_trans; - double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge - Position r_t; - if (beta_dir > 0) { - const std::array i_xyz_t {i_xyz[0]+translation_matrix(0,0), i_xyz[1]+translation_matrix(0,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-translation_matrix(0,0), i_xyz[1]-translation_matrix(0,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double beta = r_t.x * cosa + r_t.y * sina; - if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { - d = (edge - beta) / beta_dir; - if (beta_dir > 0) { - i_trans = 1; - } else { - i_trans = -1; - } - } - //gamma direction - edge = -copysign(0.5*pitch_[0], gamma_dir); - if (gamma_dir > 0) { - const std::array i_xyz_t {i_xyz[0]+translation_matrix(1,0), i_xyz[1]+translation_matrix(1,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-translation_matrix(1,0), i_xyz[1]-translation_matrix(1,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double gamma = r_t.x * cosa - r_t.y * sina; - if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { - double this_d = (edge - gamma) / gamma_dir; - if (this_d < d) { - if (gamma_dir > 0) { - i_trans = 2; - } else { - i_trans = -2; - } - d = this_d; - } - } - //y or x direction - edge = -copysign(0.5*pitch_[0], u_xy); - if (u_xy > 0) { - const std::array i_xyz_t {i_xyz[0]+translation_matrix(2,0), i_xyz[1]+translation_matrix(2,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } else { - const std::array i_xyz_t {i_xyz[0]-translation_matrix(2,0), i_xyz[1]-translation_matrix(2,1), i_xyz[2]}; - r_t = get_local_position(r, i_xyz_t); - } - double r_txy; - if (hextype==0){ - r_txy = r_t.y; - } - else { - r_txy = r_t.x; - } - if ((std::abs(r_txy - edge) > FP_PRECISION) && u_xy != 0) { - double this_d = (edge - r_txy) / u_xy; - if (this_d < d) { - if (u_xy > 0) { - i_trans = 3; - } else { - i_trans = -3; - } - d = this_d; - } - } - // Top and bottom sides - if (is_3d_) { - double z = r.z; - double z0 {copysign(0.5 * pitch_[1], u.z)}; - if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { - double this_d = (z0 - z) / u.z; - if (this_d < d) { - d = this_d; - if (u.z > 0) { - i_trans = 4; - } else { - i_trans = -4; - } - d = this_d; - } - } - } + //beta direction + double d {INFTY}; + std::array lattice_trans; + double edge = -copysign(0.5*pitch_[0], beta_dir); // Oncoming edge + double beta; + Position r_t; + if (beta_dir > 0) { + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + if (hextype==0){ + beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0;} + else{ + beta = r_t.x ;} + if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { + d = (edge - beta) / beta_dir; + if (beta_dir > 0) { + lattice_trans = {1, 0, 0}; + } else { + lattice_trans = {-1, 0, 0}; + } + } - for (int i=0;i<3;i++){ - if (i_trans < 0){ - lattice_trans[i] = -translation_matrix(-(i_trans+1),i); - }else{ - lattice_trans[i] = translation_matrix((i_trans-1),i); - } - } + // gamma direction. + edge = -copysign(0.5*pitch_[0], gamma_dir); + if (gamma_dir > 0) { + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double gamma; + if (hextype==0){ + gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0;} + else{ + gamma = r_t.x / 2.0 - r_t.y * std::sqrt(3.0) / 2.0;} + if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { + double this_d = (edge - gamma) / gamma_dir; + if (this_d < d) { + if (gamma_dir > 0) { + lattice_trans = {1, -1, 0}; + } else { + lattice_trans = {-1, 1, 0}; + } + d = this_d; + } + } - return {d, lattice_trans}; + // delta directions. + edge = -copysign(0.5*pitch_[0], delta_dir); + if (delta_dir > 0) { + const std::array i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } else { + const std::array i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; + r_t = get_local_position(r, i_xyz_t); + } + double delta; + if (hextype==0){ + delta = r_t.y ;} + else{ + delta = r_t.x / 2.0 + r_t.y * std::sqrt(3.0) / 2.0;} + if ((std::abs(delta - edge) > FP_PRECISION) && delta_dir != 0) { + double this_d = (edge - delta) / delta_dir; + if (this_d < d) { + if (delta_dir > 0) { + lattice_trans = {0, 1, 0}; + } else { + lattice_trans = {0, -1, 0}; + } + d = this_d; + } + } + + // Top and bottom sides + if (is_3d_) { + double z = r.z; + double z0 {copysign(0.5 * pitch_[1], u.z)}; + if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { + double this_d = (z0 - z) / u.z; + if (this_d < d) { + d = this_d; + if (u.z > 0) { + lattice_trans = {0, 0, 1}; + } else { + lattice_trans = {0, 0, -1}; + } + d = this_d; + } + } + } + + return {d, lattice_trans}; } //==============================================================================