diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp index 3cc8f6362d..9b64426224 100644 --- a/src/atom_vec_ellipsoid.cpp +++ b/src/atom_vec_ellipsoid.cpp @@ -536,6 +536,83 @@ void AtomVecEllipsoid::write_data_bonus(FILE *fp, int n, double *buf, int /*flag } } +/* ---------------------------------------------------------------------- + convert read_data file info from general to restricted triclinic + parent class operates on data from Velocities section of data file + child class operates on ellipsoid quaternion +------------------------------------------------------------------------- */ + +void AtomVecEllipsoid::read_data_general_to_restricted(int nlocal_previous, int nlocal) +{ + int j; + + AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal); + + // quat_g2r = quat that rotates from general to restricted triclinic + // quat_new = ellipsoid quat converted to restricted triclinic + + double quat_g2r[4],quat_new[4]; + MathExtra::mat_to_quat(domain->rotate_g2r,quat_g2r); + + for (int i = nlocal_previous; i < nlocal; i++) { + if (ellipsoid[i] < 0) continue; + j = ellipsoid[i]; + MathExtra::quatquat(quat_g2r,bonus[j].quat,quat_new); + bonus[j].quat[0] = quat_new[0]; + bonus[j].quat[1] = quat_new[1]; + bonus[j].quat[2] = quat_new[2]; + bonus[j].quat[3] = quat_new[3]; + } +} + +/* ---------------------------------------------------------------------- + convert info output by write_data from restricted to general triclinic + parent class operates on x and data from Velocities section of data file + child class operates on ellipsoid quaternion +------------------------------------------------------------------------- */ + +void AtomVecEllipsoid::write_data_restricted_to_general() +{ + AtomVec::write_data_restricted_to_general(); + + memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold"); + + for (int i = 0; i < nlocal_bonus; i++) + memcpy(quat_hold[i],bonus[i].quat,4*sizeof(double)); + + // quat_r2g = quat that rotates from restricted to general triclinic + // quat_new = ellipsoid quat converted to general triclinic + + double quat_r2g[4],quat_new[4]; + MathExtra::mat_to_quat(domain->rotate_r2g,quat_r2g); + + for (int i = 0; i < nlocal_bonus; i++) { + MathExtra::quatquat(quat_r2g,bonus[i].quat,quat_new); + bonus[i].quat[0] = quat_new[0]; + bonus[i].quat[1] = quat_new[1]; + bonus[i].quat[2] = quat_new[2]; + bonus[i].quat[3] = quat_new[3]; + } +} + +/* ---------------------------------------------------------------------- + restore info output by write_data to restricted triclinic + original data is in "hold" arrays + parent class operates on x and data from Velocities section of data file + child class operates on ellipsoid quaternion +------------------------------------------------------------------------- */ + +void AtomVecEllipsoid::write_data_restore_restricted() +{ + AtomVec::write_data_restore_restricted(); + + for (int i = 0; i < nlocal_bonus; i++) + memcpy(bonus[i].quat,quat_hold[i],4*sizeof(double)); + + memory->destroy(quat_hold); + quat_hold = nullptr; +} + /* ---------------------------------------------------------------------- set shape values in bonus data for particle I oriented aligned with xyz axes diff --git a/src/atom_vec_ellipsoid.h b/src/atom_vec_ellipsoid.h index 6e06d773fc..23c824dbf0 100644 --- a/src/atom_vec_ellipsoid.h +++ b/src/atom_vec_ellipsoid.h @@ -59,6 +59,10 @@ class AtomVecEllipsoid : public AtomVec { int pack_data_bonus(double *, int) override; void write_data_bonus(FILE *, int, double *, int) override; + void read_data_general_to_restricted(int, int); + void write_data_restricted_to_general(); + void write_data_restore_restricted(); + // unique to AtomVecEllipsoid void set_shape(int, double, double, double); @@ -69,7 +73,8 @@ class AtomVecEllipsoid : public AtomVec { int *ellipsoid; double *rmass; double **angmom; - + double **quat_hold; + int nghost_bonus, nmax_bonus; int ellipsoid_flag; double rmass_one; diff --git a/src/domain.cpp b/src/domain.cpp index 47295c338c..439dafd765 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -307,7 +307,7 @@ void Domain::set_global_box() cprime[1] = yz; cprime[2] = boxhi[2] - boxlo[2]; - // transform restricted A'B'C' to general triclinic A,B,C + // transform restricted A'B'C' to general triclinic ABC MathExtra::matvec(rotate_r2g,aprime,avec); MathExtra::matvec(rotate_r2g,bprime,bvec); diff --git a/src/math_extra.cpp b/src/math_extra.cpp index 83e548f79f..b8c9bd98df 100644 --- a/src/math_extra.cpp +++ b/src/math_extra.cpp @@ -362,6 +362,29 @@ void exyz_to_q(double *ex, double *ey, double *ez, double *q) qnormalize(q); } +/* ---------------------------------------------------------------------- + create unit quaternion from a rotation matrix + just a wrapper on exyz_to_q() + ex,ey,ez are columns of a rotation matrix +------------------------------------------------------------------------- */ + +void mat_to_quat(double mat[3][3], double *q) +{ + double ex[3],ey[3],ez[3]; + + ex[0] = mat[0][0]; + ex[1] = mat[1][0]; + ex[2] = mat[2][0]; + ey[0] = mat[0][1]; + ey[1] = mat[1][1]; + ey[2] = mat[2][1]; + ez[0] = mat[0][2]; + ez[1] = mat[1][2]; + ez[2] = mat[2][2]; + + MathExtra::exyz_to_q(ex,ey,ez,q); +} + /* ---------------------------------------------------------------------- compute space-frame ex,ey,ez from current quaternion q ex,ey,ez = space-frame coords of 1st,2nd,3rd principal axis @@ -417,6 +440,7 @@ void quat_to_mat(const double *quat, double mat[3][3]) /* ---------------------------------------------------------------------- compute rotation matrix from quaternion conjugate quat = [w i j k] + similar logic to quat_to_mat() ------------------------------------------------------------------------- */ void quat_to_mat_trans(const double *quat, double mat[3][3]) @@ -647,5 +671,4 @@ void tribbox(double *h, double radius, double *dist) /* ---------------------------------------------------------------------- */ - } diff --git a/src/math_extra.h b/src/math_extra.h index 49e128c3df..52d1d838ff 100644 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -100,6 +100,7 @@ void angmom_to_omega(double *m, double *ex, double *ey, double *ez, double *idia void omega_to_angmom(double *w, double *ex, double *ey, double *ez, double *idiag, double *m); void mq_to_omega(double *m, double *q, double *moments, double *w); void exyz_to_q(double *ex, double *ey, double *ez, double *q); +void mat_to_quat(double mat[3][3], double *quat); void q_to_exyz(double *q, double *ex, double *ey, double *ez); void quat_to_mat(const double *quat, double mat[3][3]); void quat_to_mat_trans(const double *quat, double mat[3][3]);