diff --git a/src/atom.cpp b/src/atom.cpp index 7f472e44f4..c7c13013c8 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -1190,7 +1190,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset, xdata[2] = 0.0; } - // convert atom coords from general triclinic to restricted triclinic + // convert atom coords from general to restricted triclinic // so can decide which proc owns the atom if (triclinic_general) domain->general_to_restricted_coords(xdata); diff --git a/src/atom_vec_line.cpp b/src/atom_vec_line.cpp index 95ccfe31bd..73f2e702f6 100644 --- a/src/atom_vec_line.cpp +++ b/src/atom_vec_line.cpp @@ -352,6 +352,22 @@ void AtomVecLine::data_atom_bonus(int m, const std::vector &values) double y1 = utils::numeric(FLERR, values[ivalue++], true, lmp); double x2 = utils::numeric(FLERR, values[ivalue++], true, lmp); double y2 = utils::numeric(FLERR, values[ivalue++], true, lmp); + + // convert x1/y1 and x2/y2 from general to restricted triclniic + // x is already restricted triclinic + + if (domain->triclinic_general) { + double coords[3]; + coords[0] = x1; coords[1] = y1; coords[2] = 0.0; + domain->general_to_restricted_coords(coords); + x1 = coords[0]; y1 = coords[1]; + coords[0] = x2; coords[1] = y2; coords[2] = 0.0; + domain->general_to_restricted_coords(coords); + x2 = coords[0]; y2 = coords[1]; + } + + // calculate length and theta + double dx = x2 - x1; double dy = y2 - y1; double length = sqrt(dx * dx + dy * dy); @@ -437,29 +453,6 @@ void AtomVecLine::data_atom_post(int ilocal) omega[ilocal][2] = 0.0; } -/* ---------------------------------------------------------------------- - 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 bonus theta -------------------------------------------------------------------------- */ - -void AtomVecLine::read_data_general_to_restricted(int nlocal_previous, int nlocal) -{ - AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal); - - double btheta; - double theta_g2r = 2.0*acos(domain->quat_g2r[0]); - - for (int i = nlocal_previous; i < nlocal; i++) { - if (line[i] < 0) continue; - btheta = bonus[line[i]].theta; - btheta += theta_g2r; - if (btheta > MathConst::MY_PI) btheta -= MathConst::MY_2PI; - else if (btheta <= -MathConst::MY_PI) btheta += MathConst::MY_2PI; - bonus[line[i]].theta = btheta; - } -} - /* ---------------------------------------------------------------------- modify values for AtomVec::pack_data() to pack ------------------------------------------------------------------------- */ @@ -501,8 +494,14 @@ int AtomVecLine::pack_data_bonus(double *buf, int /*flag*/) int i, j; double length, theta; double xc, yc, x1, x2, y1, y2; + double coords[3]; + + int triclinic_general = domain->triclinic_general; + + double **x_bonus; + if (triclinic_general) x_bonus = x_hold; + else x_bonus = x; - double **x = atom->x; tagint *tag = atom->tag; int nlocal = atom->nlocal; @@ -514,8 +513,9 @@ int AtomVecLine::pack_data_bonus(double *buf, int /*flag*/) j = line[i]; length = bonus[j].length; theta = bonus[j].theta; - xc = x[i][0]; - yc = x[i][1]; + + xc = x_bonus[i][0]; + yc = x_bonus[i][1]; x1 = xc - 0.5 * cos(theta) * length; y1 = yc - 0.5 * sin(theta) * length; x2 = xc + 0.5 * cos(theta) * length; @@ -524,6 +524,20 @@ int AtomVecLine::pack_data_bonus(double *buf, int /*flag*/) buf[m++] = y1; buf[m++] = x2; buf[m++] = y2; + + // if triclinic_general: + // rotate 4 buf values from restricted to general triclinic + // output by write_data_bonus() as x1/y1 and x2/y2 + + if (triclinic_general) { + coords[0] = buf[m-4]; coords[1] = buf[m-3]; coords[2] = 0.0; + domain->restricted_to_general_coords(coords); + buf[m-4] = coords[0]; buf[m-3] = coords[1]; + coords[0] = buf[m-2]; coords[1] = buf[m-1]; coords[2] = 0.0; + domain->restricted_to_general_coords(coords); + buf[m-2] = coords[0]; buf[m-1] = coords[1]; + } + } else m += size_data_bonus; } diff --git a/src/atom_vec_line.h b/src/atom_vec_line.h index 2503f55d8d..740c541916 100644 --- a/src/atom_vec_line.h +++ b/src/atom_vec_line.h @@ -53,7 +53,6 @@ class AtomVecLine : public AtomVec { void create_atom_post(int) override; void data_atom_post(int) override; - void read_data_general_to_restricted(int, int); void pack_data_pre(int) override; void pack_data_post(int) override; diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp index f83a282c35..b0849e6717 100644 --- a/src/atom_vec_tri.cpp +++ b/src/atom_vec_tri.cpp @@ -510,8 +510,18 @@ void AtomVecTri::data_atom_bonus(int m, const std::vector &values) MathExtra::sub3(c3, c1, c3mc1); double size = MAX(MathExtra::len3(c2mc1), MathExtra::len3(c3mc1)); + // convert c1,c2,c3 from general to restricted triclniic + // x is already restricted triclinic + + if (domain->triclinic_general) { + domain->general_to_restricted_coords(c1); + domain->general_to_restricted_coords(c2); + domain->general_to_restricted_coords(c3); + } + // centroid = 1/3 of sum of vertices - + // error if centroid is not within EPSILON of Atoms section coord + double centroid[3]; centroid[0] = (c1[0] + c2[0] + c3[0]) / 3.0; centroid[1] = (c1[1] + c2[1] + c3[1]) / 3.0; @@ -659,90 +669,6 @@ void AtomVecTri::data_atom_post(int ilocal) angmom[ilocal][2] = 0.0; } -/* ---------------------------------------------------------------------- - 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 bonus quat -------------------------------------------------------------------------- */ - -void AtomVecTri::read_data_general_to_restricted(int nlocal_previous, int nlocal) -{ - AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal); - - double quat[4]; - double *bquat; - - for (int i = nlocal_previous; i < nlocal; i++) { - if (tri[i] < 0) continue; - bquat = bonus[tri[i]].quat; - MathExtra::quatquat(domain->quat_g2r,bquat,quat); - bquat[0] = quat[0]; - bquat[1] = quat[1]; - bquat[2] = quat[2]; - bquat[3] = quat[3]; - MathExtra::qnormalize(bquat); - } -} - -/* ---------------------------------------------------------------------- - 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 bonus quat -------------------------------------------------------------------------- */ - -void AtomVecTri::write_data_restricted_to_general() -{ - AtomVec::write_data_restricted_to_general(); - - double quat[4]; - double *bquat; - double *quat_r2g = domain->quat_r2g; - - memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold"); - - for (int i = 0; i < nlocal_bonus; i++) { - bquat = bonus[i].quat; - quat_hold[i][0] = bquat[0]; - quat_hold[i][1] = bquat[1]; - quat_hold[i][2] = bquat[2]; - quat_hold[i][3] = bquat[3]; - - MathExtra::quatquat(quat_r2g,bquat,quat); - bquat[0] = quat[0]; - bquat[1] = quat[1]; - bquat[2] = quat[2]; - bquat[3] = quat[3]; - MathExtra::qnormalize(bquat); - } -} - -/* ---------------------------------------------------------------------- - 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 bonus quat -------------------------------------------------------------------------- */ - -void AtomVecTri::write_data_restore_restricted() -{ - AtomVec::write_data_restore_restricted(); - - if (!quat_hold) return; - - double *bquat; - - for (int i = 0; i < nlocal_bonus; i++) { - bquat = bonus[i].quat; - bquat[0] = quat_hold[i][0]; - bquat[1] = quat_hold[i][1]; - bquat[2] = quat_hold[i][2]; - bquat[3] = quat_hold[i][3]; - } - - memory->destroy(quat_hold); - quat_hold = nullptr; -} - /* ---------------------------------------------------------------------- modify values for AtomVec::pack_data() to pack ------------------------------------------------------------------------- */ @@ -792,7 +718,12 @@ int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/) double dc1[3], dc2[3], dc3[3]; double p[3][3]; - double **x = atom->x; + int triclinic_general = domain->triclinic_general; + + double **x_bonus; + if (triclinic_general) x_bonus = x_hold; + else x_bonus = x; + tagint *tag = atom->tag; int nlocal = atom->nlocal; @@ -806,9 +737,10 @@ int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/) MathExtra::matvec(p, bonus[j].c1, dc1); MathExtra::matvec(p, bonus[j].c2, dc2); MathExtra::matvec(p, bonus[j].c3, dc3); - xc = x[i][0]; - yc = x[i][1]; - zc = x[i][2]; + + xc = x_bonus[i][0]; + yc = x_bonus[i][1]; + zc = x_bonus[i][2]; buf[m++] = xc + dc1[0]; buf[m++] = yc + dc1[1]; buf[m++] = zc + dc1[2]; @@ -818,6 +750,17 @@ int AtomVecTri::pack_data_bonus(double *buf, int /*flag*/) buf[m++] = xc + dc3[0]; buf[m++] = yc + dc3[1]; buf[m++] = zc + dc3[2]; + + // if triclinic_general: + // rotate 9 buf values from restricted to general triclinic + // output by write_data_bonus() as c1,c2,c3 + + if (triclinic_general) { + domain->restricted_to_general_coords(&buf[m-9]); + domain->restricted_to_general_coords(&buf[m-6]); + domain->restricted_to_general_coords(&buf[m-3]); + } + } else m += size_data_bonus; } diff --git a/src/atom_vec_tri.h b/src/atom_vec_tri.h index 91a7a20e6d..5a3b831c0d 100644 --- a/src/atom_vec_tri.h +++ b/src/atom_vec_tri.h @@ -55,9 +55,6 @@ class AtomVecTri : public AtomVec { void create_atom_post(int) override; void data_atom_post(int) override; - void read_data_general_to_restricted(int, int); - void write_data_restricted_to_general(); - void write_data_restore_restricted(); void pack_data_pre(int) override; void pack_data_post(int) override; diff --git a/src/lmprestart.h b/src/lmprestart.h index 43ef5766ac..73eb5e618e 100644 --- a/src/lmprestart.h +++ b/src/lmprestart.h @@ -38,7 +38,7 @@ enum{VERSION,SMALLINT,TAGINT,BIGINT, EXTRA_BOND_PER_ATOM,EXTRA_ANGLE_PER_ATOM,EXTRA_DIHEDRAL_PER_ATOM, EXTRA_IMPROPER_PER_ATOM,EXTRA_SPECIAL_PER_ATOM,ATOM_MAXSPECIAL, NELLIPSOIDS,NLINES,NTRIS,NBODIES,ATIME,ATIMESTEP,LABELMAP, - TRICLINIC_GENERAL,TRICLINIC_GENERAL_FLIP,QUAT_G2R}; + TRICLINIC_GENERAL,TRICLINIC_GENERAL_FLIP,ROTATE_G2R}; #define LB_FACTOR 1.1 diff --git a/src/read_restart.cpp b/src/read_restart.cpp index ca5eb25762..2224b9eeac 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -137,21 +137,6 @@ void ReadRestart::command(int narg, char **arg) atom->avec->grow(n); n = atom->nmax; - // setup simulation box - // for general triclinic, need to generate additional info which - // Domain::define_general_triclinic() would have created - - if (domain->triclinic_general) { - MathExtra::qconjugate(domain->quat_g2r,domain->quat_r2g); - MathExtra::quat_to_mat(domain->quat_g2r,domain->rotate_g2r); - if (domain->triclinic_general_flip) { - domain->rotate_g2r[2][0] = -domain->rotate_g2r[2][0]; - domain->rotate_g2r[2][1] = -domain->rotate_g2r[2][1]; - domain->rotate_g2r[2][2] = -domain->rotate_g2r[2][2]; - } - MathExtra::transpose3(domain->rotate_g2r,domain->rotate_r2g); - } - domain->print_box(" "); domain->set_initial_box(0); domain->set_global_box(); @@ -803,8 +788,9 @@ void ReadRestart::header() domain->triclinic_general = read_int(); } else if (flag == TRICLINIC_GENERAL_FLIP) { domain->triclinic_general_flip = read_int(); - } else if (flag == QUAT_G2R) { - read_double_vec(4,domain->quat_g2r); + } else if (flag == ROTATE_G2R) { + read_double_vec(9,&domain->rotate_g2r[0][0]); + MathExtra::transpose3(domain->rotate_g2r,domain->rotate_r2g); } else if (flag == SPECIAL_LJ) { read_int(); diff --git a/src/write_restart.cpp b/src/write_restart.cpp index 16a59ebc20..e265ceadc1 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -451,7 +451,7 @@ void WriteRestart::header() write_int(TRICLINIC_GENERAL,domain->triclinic_general); if (domain->triclinic_general) { write_int(TRICLINIC_GENERAL_FLIP,domain->triclinic_general_flip); - write_double_vec(QUAT_G2R,4,domain->quat_g2r); + write_double_vec(ROTATE_G2R,9,&domain->rotate_g2r[0][0]); } write_double_vec(SPECIAL_LJ,3,&force->special_lj[1]);