more debugging for fcc example

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);

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);

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);