diff --git a/src/create_box.cpp b/src/create_box.cpp index 99dd74d566..e42c7c32f0 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -96,7 +96,7 @@ void CreateBox::command(int narg, char **arg) // setup general triclinic box (with no region) // read next box extent arguments to create ABC edge vectors + origin - // set_general_triclinic() converts + // setup_general_triclinic() converts // ABC edge vectors + origin to restricted triclinic } else if (triclinic_general) { @@ -117,7 +117,7 @@ void CreateBox::command(int narg, char **arg) // use lattice2box() to generate origin and ABC vectors // origin = abc lo - // ABC vector = hi in one dim - orign + // ABC vectors = hi in one dim - orign double avec[3],bvec[3],cvec[3],origin[3]; double px,py,pz; @@ -129,12 +129,10 @@ void CreateBox::command(int narg, char **arg) origin[2] = pz; px = ahi; py = blo; pz = clo; - printf("CB PXYZ %g %g %g\n",px,py,pz); domain->lattice->lattice2box(px,py,pz); avec[0] = px - origin[0]; avec[1] = py - origin[1]; avec[2] = pz - origin[2]; - printf("CB AVEC %g %g %g\n",avec[0],avec[1],avec[2]); px = alo; py = bhi; pz = clo; domain->lattice->lattice2box(px,py,pz); @@ -150,7 +148,7 @@ void CreateBox::command(int narg, char **arg) // setup general triclinic box within Domain - domain->set_general_triclinic(avec,bvec,cvec,origin); + domain->setup_general_triclinic(avec,bvec,cvec,origin); } // if molecular, zero out topology info diff --git a/src/domain.cpp b/src/domain.cpp index f1b0c5a20e..29e0e2cd89 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -527,8 +527,8 @@ void Domain::reset_box() create rotation matrices for general <--> restricted transformations ------------------------------------------------------------------------- */ -void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, - double *cvec_caller, double *origin_caller) +void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller, + double *cvec_caller, double *origin_caller) { if (triclinic || triclinic_general) error->all(FLERR,"General triclinic box edge vectors are already set"); @@ -547,15 +547,15 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, cvec[1] = cvec_caller[1]; cvec[2] = cvec_caller[2]; - tri_origin[0] = origin_caller[0]; - tri_origin[1] = origin_caller[1]; - tri_origin[2] = origin_caller[2]; + gtri_origin[0] = origin_caller[0]; + gtri_origin[1] = origin_caller[1]; + gtri_origin[2] = origin_caller[2]; // error check on cvec for 2d systems if (dimension == 2 && (cvec[0] != 0.0 || cvec[1] != 0.0)) error->all(FLERR,"General triclinic box edge vector C invalid for 2d system"); - + // error check for co-planar A,B,C double abcross[3]; @@ -563,8 +563,8 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, double dot = MathExtra::dot3(abcross,cvec); if (dot == 0.0) error->all(FLERR,"General triclinic box edge vectors are co-planar"); - - // quat1 = convert A into A' along x-axis + + // quat1 = convert A into A' along +x-axis // rot1 = unit vector to rotate A around // theta1 = angle of rotation calculated from // A dot xunit = Ax = |A| cos(theta1) @@ -583,45 +583,34 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, double rotmat1[3][3]; MathExtra::quat_to_mat(quat1,rotmat1); - // DEBUG - double afinal[3]; - MathExtra::matvec(rotmat1,avec,afinal); - printf("AFINAL %g %g %g\n",afinal[0],afinal[1],afinal[2]); - // B1 = rotation of B by quat1 rotation matrix double bvec1[3]; MathExtra::matvec(rotmat1,bvec,bvec1); - printf("BVEC1 %g %g %g\n",bvec1[0],bvec1[1],bvec1[2]); - // quat2 = rotation to convert B1 into B' in xy plane // Byz1 = projection of B1 into yz plane - // rot2 = unit vector to rotate B1 around = -x axis + // +xaxis = unit vector to rotate B1 around // theta2 = angle of rotation calculated from // Byz1 dot yunit = B1y = |Byz1| cos(theta2) - + // theta2 via acos() is positive (0 to PI) + // positive is valid if B1z < 0.0 else flip sign of theta2 + double byzvec1[3],quat2[4]; MathExtra::copy3(bvec1,byzvec1); byzvec1[0] = 0.0; double byzvec1_len = MathExtra::len3(byzvec1); - double rot2[3] = {-1.0, 0.0, 0.0}; - if (byzvec1[2] < 0.0) rot2[0] = 1.0; double theta2 = acos(bvec1[1]/byzvec1_len); - MathExtra::axisangle_to_quat(rot2,theta2,quat2); + if (bvec1[2] > 0.0) theta2 = -theta2; + MathExtra::axisangle_to_quat(xaxis,theta2,quat2); - // DEBUG - double rotmat2[3][3]; - MathExtra::quat_to_mat(quat2,rotmat2); - double bfinal[3]; - MathExtra::matvec(rotmat1,bvec1,bfinal); - printf("BFINAL %g %g %g\n",bfinal[0],bfinal[1],bfinal[2]); - - // quat = product of quat2 * quat1 = transformation via single quat + // quat = transformation via single quat = quat2 * quat1 // rotate_g2r = general to restricted transformation matrix + // if dot < 0.0 (A x B not in C direction) + // flip sign of z component of transform, + // by flipping sign of 3rd row of rotate_g2r matrix // rotate_r2g = restricted to general transformation matrix - // if A x B not in direction of C, flip sign of z component of transform - // done by flipping sign of 3rd row of rotate_g2r matrix + // simply a transpose of rotate_g2r since orthonormal double quat[4]; MathExtra::quatquat(quat2,quat1,quat); @@ -644,9 +633,9 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, // set restricted triclinic boxlo, boxhi, and tilt factors - boxlo[0] = tri_origin[0]; - boxlo[1] = tri_origin[1]; - boxlo[2] = tri_origin[2]; + boxlo[0] = gtri_origin[0]; + boxlo[1] = gtri_origin[1]; + boxlo[2] = gtri_origin[2]; boxhi[0] = boxlo[0] + aprime[0]; boxhi[1] = boxlo[1] + bprime[1]; @@ -658,14 +647,10 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, // debug - printf("Cvec: %g %g %g\n",cvec[0],cvec[1],cvec[2]); - printf("ABcross: %g %g %g\n",abcross[0],abcross[1],abcross[2]); - printf("Dot: %g\n",dot); - printf("Avec: %g %g %g\n",avec[0],avec[1],avec[2]); printf("Theta1: %g\n",theta1); printf("Rotvec1: %g %g %g\n",rot1[0],rot1[1],rot1[2]); printf("Theta2: %g\n",theta2); - printf("Rotvec2: %g %g %g\n",rot2[0],rot2[1],rot2[2]); + printf("Rotvec2: %g %g %g\n",xaxis[0],xaxis[1],xaxis[2]); printf("Quat: %g %g %g %g\n",quat[0],quat[1],quat[2],quat[3]); double angle = 2.0*acos(quat[0]); printf("Theta: %g\n",angle); @@ -676,14 +661,6 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, printf("Length A: %g %g\n",MathExtra::len3(avec),MathExtra::len3(aprime)); printf("Length B: %g %g\n",MathExtra::len3(bvec),MathExtra::len3(bprime)); printf("Length C: %g %g\n",MathExtra::len3(cvec),MathExtra::len3(cprime)); - - double coord1[3] = {0.5,0.0,0.0}; - double coord2[3] = {0.5,0.0,0.3}; - double newcoord[3]; - MathExtra::matvec(rotate_g2r,coord1,newcoord); - printf("Atom1: %g %g %g\n",newcoord[0],newcoord[1],newcoord[2]); - MathExtra::matvec(rotate_g2r,coord2,newcoord); - printf("Atom2: %g %g %g\n",newcoord[0],newcoord[1],newcoord[2]); } /* ---------------------------------------------------------------------- @@ -695,9 +672,9 @@ void Domain::general_to_restricted(double *x) double xnew[3]; MathExtra::matvec(rotate_g2r,x,xnew); - x[0] = xnew[0] + tri_origin[0]; - x[1] = xnew[1] + tri_origin[1]; - x[2] = xnew[2] + tri_origin[2]; + x[0] = xnew[0] + gtri_origin[0]; + x[1] = xnew[1] + gtri_origin[1]; + x[2] = xnew[2] + gtri_origin[2]; } /* ---------------------------------------------------------------------- @@ -708,9 +685,9 @@ void Domain::restricted_to_general(double *x) { double xshift[3],xnew[3]; - xshift[0] = x[0] - tri_origin[0]; - xshift[1] = x[1] - tri_origin[1]; - xshift[2] = x[2] - tri_origin[2]; + xshift[0] = x[0] - gtri_origin[0]; + xshift[1] = x[1] - gtri_origin[1]; + xshift[2] = x[2] - gtri_origin[2]; MathExtra::matvec(rotate_r2g,xshift,xnew); x[0] = xnew[0]; x[1] = xnew[1]; diff --git a/src/domain.h b/src/domain.h index e77af5d968..d09f4cf70d 100644 --- a/src/domain.h +++ b/src/domain.h @@ -91,7 +91,7 @@ class Domain : protected Pointers { // general triclinic box double avec[3], bvec[3], cvec[3]; // ABC edge vectors of general triclinic box - double tri_origin[3]; // origin of general triclinic box + double gtri_origin[3]; // origin of general triclinic box double rotate_g2r[3][3]; // rotation matrix from general --> restricted tri double rotate_r2g[3][3]; // rotation matrix from restricted --> general tri @@ -141,7 +141,7 @@ class Domain : protected Pointers { void image_flip(int, int, int); int ownatom(int, double *, imageint *, int); - void set_general_triclinic(double *, double *, double *, double *); + void setup_general_triclinic(double *, double *, double *, double *); void general_to_restricted(double *); void restricted_to_general(double *); diff --git a/src/read_data.cpp b/src/read_data.cpp index bd998ed3f5..75d1c23b75 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -480,7 +480,7 @@ void ReadData::command(int narg, char **arg) ellipsoidflag = lineflag = triflag = bodyflag = 0; xloxhi_flag = yloyhi_flag = zlozhi_flag = tilt_flag = 0; - avec_flag = bvec_flag = cvec_flag = 0; + avec_flag = bvec_flag = cvec_flag = gtri_origin_flag = 0; // values in this data file @@ -495,7 +495,7 @@ void ReadData::command(int narg, char **arg) bvec[0] = bvec[1] = bvec[2] = 0.0; cvec[0] = cvec[1] = cvec[2] = 0.0; avec[0] = bvec[1] = cvec[2] = 1.0; - tri_origin[0] = tri_origin[1] = tri_origin[2] = 0.0; + gtri_origin[0] = gtri_origin[1] = gtri_origin[2] = 0.0; keyword[0] = '\0'; @@ -518,7 +518,7 @@ void ReadData::command(int narg, char **arg) // check if simulation box specified consistently - if (!avec_flag && !bvec_flag && !cvec_flag) { + if (!avec_flag && !bvec_flag && !cvec_flag && !gtri_origin_flag) { triclinic = triclinic_general = 0; if (tilt_flag) triclinic = 1; } else { @@ -579,11 +579,11 @@ void ReadData::command(int narg, char **arg) } // general triclinic box - // set_general_triclinic() converts - // ABC edge vectors + origin to restricted triclinic + // setup_general_triclinic() converts + // ABC edge vectors + gtri_origin to restricted triclinic } else if (triclinic_general) { - domain->set_general_triclinic(avec,bvec,cvec,tri_origin); + domain->setup_general_triclinic(avec,bvec,cvec,gtri_origin); } } @@ -609,8 +609,8 @@ void ReadData::command(int narg, char **arg) errflag = 1; if (cvec[0] != domain->cvec[0] || cvec[1] != domain->cvec[1] || cvec[2] != domain->cvec[2]) errflag = 1; - if (tri_origin[0] != domain->tri_origin[0] || tri_origin[1] != domain->tri_origin[1] || - tri_origin[2] != domain->tri_origin[2]) + if (gtri_origin[0] != domain->gtri_origin[0] || gtri_origin[1] != domain->gtri_origin[1] || + gtri_origin[2] != domain->gtri_origin[2]) errflag = 1; if (errflag) error->all(FLERR,"Read_data subsequent file ABC vectors must be same as first file"); @@ -1414,26 +1414,29 @@ void ReadData::header(int firstpass) xz = utils::numeric(FLERR, words[1], false, lmp); yz = utils::numeric(FLERR, words[2], false, lmp); - } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\f+\\s+\\avec\\s")) { + } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\avec\\s")) { avec_flag = 1; avec[0] = utils::numeric(FLERR, words[0], false, lmp); avec[1] = utils::numeric(FLERR, words[1], false, lmp); avec[2] = utils::numeric(FLERR, words[2], false, lmp); - tri_origin[0] = utils::numeric(FLERR, words[3], false, lmp); - } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\f+\\s+\\bvec\\s")) { + } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\bvec\\s")) { bvec_flag = 1; bvec[0] = utils::numeric(FLERR, words[0], false, lmp); bvec[1] = utils::numeric(FLERR, words[1], false, lmp); bvec[2] = utils::numeric(FLERR, words[2], false, lmp); - tri_origin[1] = utils::numeric(FLERR, words[3], false, lmp); - } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\f+\\s+\\cvec\\s")) { + } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\cvec\\s")) { cvec_flag = 1; cvec[0] = utils::numeric(FLERR, words[0], false, lmp); cvec[1] = utils::numeric(FLERR, words[1], false, lmp); cvec[2] = utils::numeric(FLERR, words[2], false, lmp); - tri_origin[2] = utils::numeric(FLERR, words[3], false, lmp); + + } else if (utils::strmatch(line, "^\\s*\\f+\\s+\\f+\\s+\\f+\\s+\\gtri\\s+origin\\s")) { + gtri_origin_flag = 1; + gtri_origin[0] = utils::numeric(FLERR, words[0], false, lmp); + gtri_origin[1] = utils::numeric(FLERR, words[1], false, lmp); + gtri_origin[2] = utils::numeric(FLERR, words[2], false, lmp); } else break; diff --git a/src/read_data.h b/src/read_data.h index 6439fbd096..e3c269e685 100644 --- a/src/read_data.h +++ b/src/read_data.h @@ -63,10 +63,10 @@ class ReadData : public Command { double boxlo[3], boxhi[3]; double xy, xz, yz; - double avec[3], bvec[3], cvec[3], tri_origin[3]; + double avec[3], bvec[3], cvec[3], gtri_origin[3]; int triclinic, triclinic_general; int xloxhi_flag, yloyhi_flag, zlozhi_flag, tilt_flag; - int avec_flag, bvec_flag, cvec_flag; + int avec_flag, bvec_flag, cvec_flag, gtri_origin_flag; // optional args diff --git a/src/write_data.cpp b/src/write_data.cpp index 03e09cf462..ad0c508d15 100644 --- a/src/write_data.cpp +++ b/src/write_data.cpp @@ -331,10 +331,12 @@ void WriteData::header() fmt::print(fp,"{} {} {} xy xz yz\n",domain->xy,domain->xz,domain->yz); } else if (triclinic_general) { - fmt::print(fp,"\n{} {} {} {} avec\n{} {} {} {} bvec\n{} {} {} {} cvec\n", - domain->avec[0],domain->avec[1],domain->avec[2],domain->tri_origin[0], - domain->bvec[0],domain->bvec[1],domain->bvec[2],domain->tri_origin[1], - domain->cvec[0],domain->cvec[1],domain->cvec[2],domain->tri_origin[2]); + fmt::print(fp,"\n{} {} {} avec\n{} {} {} bvec\n{} {} {} cvec\n", + domain->avec[0],domain->avec[1],domain->avec[2], + domain->bvec[0],domain->bvec[1],domain->bvec[2], + domain->cvec[0],domain->cvec[1],domain->cvec[2]); + fmt::print(fp,"{} {} {} gtri origin\n", + domain->gtri_origin[0],domain->gtri_origin[1],domain->gtri_origin[2]); } }