ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

code cleanup

Browse Source
This commit is contained in:
Steve Plimpton
2023-09-02 18:45:37 -06:00
parent 8c3ab47fd6
commit 57f6526e53
6 changed files with 59 additions and 79 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/create_box.cpp
View File

@ -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

73
src/domain.cpp
View File

@ -527,7 +527,7 @@ void Domain::reset_box()
create rotation matrices for general <--> restricted transformations
------------------------------------------------------------------------- */
void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller,
void Domain::setup_general_triclinic(double *avec_caller, double *bvec_caller,
double *cvec_caller, double *origin_caller)
{
if (triclinic || triclinic_general)
@ -547,9 +547,9 @@ 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
@ -564,7 +564,7 @@ void Domain::set_general_triclinic(double *avec_caller, double *bvec_caller,
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];

4
src/domain.h
View File

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

31
src/read_data.cpp
View File

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

4
src/read_data.h
View File

@ -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

10
src/write_data.cpp
View File

@ -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]);
}
}