diff --git a/src/create_box.cpp b/src/create_box.cpp index 24dff6c01f..99dd74d566 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -129,10 +129,12 @@ 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); diff --git a/src/domain.cpp b/src/domain.cpp index 69ca63169c..f1b0c5a20e 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -577,32 +577,46 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller, MathExtra::norm3(rot1); double theta1 = acos(avec[0]/avec_len); MathExtra::axisangle_to_quat(rot1,theta1,quat1); - + // rotmat1 = rotation matrix associated with quat1 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 // theta2 = angle of rotation calculated from - // Byz1 dot yunit = B1y = |Byz1} cos(theta2) + // Byz1 dot yunit = B1y = |Byz1| cos(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); + // 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 // rotate_g2r = general to restricted transformation matrix // rotate_r2g = restricted to general transformation matrix @@ -644,6 +658,14 @@ 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("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); diff --git a/src/lattice.cpp b/src/lattice.cpp index 2e7e740fb8..a5dc9714cf 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -186,6 +186,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) a1[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); a1[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); iarg += 4; + } else if (strcmp(arg[iarg],"a2") == 0) { if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "lattice a2", error); if (style != CUSTOM) @@ -254,11 +255,12 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) // reset scale for LJ units (input scale is rho*) // scale = (Nbasis/(Vprimitive * rho*)) ^ (1/dim) - + // use fabs() in case a1,a2,a3 are not right-handed for general triclinic + if (strcmp(update->unit_style,"lj") == 0) { double vec[3]; cross(a2,a3,vec); - double volume = dot(a1,vec); + double volume = fabs(dot(a1,vec)); scale = pow(nbasis/volume/scale,1.0/dimension); } @@ -497,7 +499,7 @@ void Lattice::setup_transform() ------------------------------------------------------------------------- */ void Lattice::lattice2box(double &x, double &y, double &z) -{ +{ double x1 = primitive[0][0]*x + primitive[0][1]*y + primitive[0][2]*z; double y1 = primitive[1][0]*x + primitive[1][1]*y + primitive[1][2]*z; double z1 = primitive[2][0]*x + primitive[2][1]*y + primitive[2][2]*z; @@ -590,7 +592,7 @@ void Lattice::cross(double *x, double *y, double *z) void Lattice::bbox(int flag, double x, double y, double z, double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) -{ +{ if (flag == 0) lattice2box(x,y,z); else box2lattice(x,y,z);