From c5b2d662833a80a5bf2fb8b8294292313872f475 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 15 Sep 2023 16:49:56 -0600 Subject: [PATCH] upgrades to read/write data commands --- src/atom_vec_body.cpp | 6 +++--- src/atom_vec_ellipsoid.cpp | 6 +++--- src/atom_vec_tri.cpp | 4 ++-- src/domain.cpp | 24 +++++++++++++++--------- src/domain.h | 5 +++-- src/dump_atom.cpp | 5 +++-- src/lmprestart.h | 2 +- src/read_restart.cpp | 24 ++++++++++++++++++++---- src/write_restart.cpp | 4 ++-- 9 files changed, 52 insertions(+), 28 deletions(-) diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp index 7721e1540b..898b065749 100644 --- a/src/atom_vec_body.cpp +++ b/src/atom_vec_body.cpp @@ -586,11 +586,11 @@ void AtomVecBody::write_data_restricted_to_general() { AtomVec::write_data_restricted_to_general(); - double quat[4],quat_r2g[4]; + double quat[4]; double *bquat; - + double *quat_r2g = domain->quat_r2g; + memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold"); - MathExtra::qconjugate(domain->quat_g2r,quat_r2g); for (int i = 0; i < nlocal_bonus; i++) { bquat = bonus[i].quat; diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp index ceb813cbe2..d51f7af78f 100644 --- a/src/atom_vec_ellipsoid.cpp +++ b/src/atom_vec_ellipsoid.cpp @@ -495,11 +495,11 @@ void AtomVecEllipsoid::write_data_restricted_to_general() { AtomVec::write_data_restricted_to_general(); - double quat[4],quat_r2g[4]; + double quat[4]; double *bquat; - + double *quat_r2g = domain->quat_r2g; + memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold"); - MathExtra::qconjugate(domain->quat_g2r,quat_r2g); for (int i = 0; i < nlocal_bonus; i++) { bquat = bonus[i].quat; diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp index 7dfb50093b..f83a282c35 100644 --- a/src/atom_vec_tri.cpp +++ b/src/atom_vec_tri.cpp @@ -694,11 +694,11 @@ void AtomVecTri::write_data_restricted_to_general() { AtomVec::write_data_restricted_to_general(); - double quat[4],quat_r2g[4]; + double quat[4]; double *bquat; + double *quat_r2g = domain->quat_r2g; memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold"); - MathExtra::qconjugate(domain->quat_g2r,quat_r2g); for (int i = 0; i < nlocal_bonus; i++) { bquat = bonus[i].quat; diff --git a/src/domain.cpp b/src/domain.cpp index e6ef232a6d..68eafab356 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -580,13 +580,18 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller, error->all(FLERR,"General triclinic box edge vector C invalid for 2d system"); // error check for co-planar A,B,C + // if A x B is in C direction, only a rotation to restricted tri is needed + // since restricted tri box will obey right-hand rule + // if not, an inversion (flip) of C vector is also needed double abcross[3]; MathExtra::cross3(avec,bvec,abcross); double dot = MathExtra::dot3(abcross,cvec); if (dot == 0.0) error->all(FLERR,"General triclinic box edge vectors are co-planar"); - + if (dot > 0.0) triclinic_general_flip = 0; + else triclinic_general_flip = 1; + // quat1 = convert A into A' along +x-axis // rot1 = unit vector to rotate A around // theta1 = angle of rotation calculated from @@ -627,18 +632,19 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller, if (bvec1[2] > 0.0) theta2 = -theta2; MathExtra::axisangle_to_quat(xaxis,theta2,quat2); - // quat_g2r = 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 + // quat_g2r = rotation via single quat = quat2 * quat1 + // quat_r2g = rotation from restricted to general + // rotate_g2r = general to restricted rotation matrix + // include flip of C vector if needed to obey right-hand rule + // rotate_r2g = restricted to general rotation matrix // simply a transpose of rotate_g2r since orthonormal MathExtra::quatquat(quat2,quat1,quat_g2r); + MathExtra::qconjugate(quat_g2r,quat_r2g); + MathExtra::quat_to_mat(quat_g2r,rotate_g2r); - if (dot < 0.0) { + if (triclinic_general_flip) { rotate_g2r[2][0] = -rotate_g2r[2][0]; rotate_g2r[2][1] = -rotate_g2r[2][1]; rotate_g2r[2][2] = -rotate_g2r[2][2]; @@ -646,7 +652,7 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller, MathExtra::transpose3(rotate_g2r,rotate_r2g); - // transform general ABC to restricted triclinic A'B'C' + // rotate general ABC to restricted triclinic A'B'C' double aprime[3],bprime[3],cprime[3]; MathExtra::matvec(rotate_g2r,avec,aprime); diff --git a/src/domain.h b/src/domain.h index 7c5510d4c2..61a88724e7 100644 --- a/src/domain.h +++ b/src/domain.h @@ -41,6 +41,7 @@ class Domain : protected Pointers { int triclinic; // 0 = orthog box, 1 = triclinic (restricted or general) int triclinic_general; // 1 if mapping to/from general triclinic is stored, 0 if not + int triclinic_general_flip; // 1 if general tri rotation needs to invert C edge vector // orthogonal box @@ -90,9 +91,9 @@ class Domain : protected Pointers { // general triclinic box // boxlo = lower left corner - + double avec[3], bvec[3], cvec[3]; // ABC edge vectors of general triclinic box - double quat_g2r[4]; // quaternion for general --> restricted rotation + double quat_g2r[4], quat_r2g[4]; // quaternions for general <--> restricted rotations double rotate_g2r[3][3]; // rotation matrix from general --> restricted tri double rotate_r2g[3][3]; // rotation matrix from restricted --> general tri diff --git a/src/dump_atom.cpp b/src/dump_atom.cpp index dd46452d9e..151446bdfd 100644 --- a/src/dump_atom.cpp +++ b/src/dump_atom.cpp @@ -334,13 +334,14 @@ void DumpAtom::header_binary_triclinic_general(bigint ndump) fwrite(&update->ntimestep,sizeof(bigint),1,fp); fwrite(&ndump,sizeof(bigint),1,fp); - int general_tri = 2; - fwrite(&general_tri,sizeof(int),1,fp); + int triclinic_general_flag = 2; + fwrite(&triclinic_general_flag,sizeof(int),1,fp); fwrite(&domain->boundary[0][0],6*sizeof(int),1,fp); fwrite(domain->avec,3*sizeof(double),1,fp); fwrite(domain->bvec,3*sizeof(double),1,fp); fwrite(domain->cvec,3*sizeof(double),1,fp); fwrite(domain->boxlo,3*sizeof(double),1,fp); + fwrite(&size_one,sizeof(int),1,fp); header_unit_style_binary(); header_time_binary(); diff --git a/src/lmprestart.h b/src/lmprestart.h index 2ed1d7db11..43ef5766ac 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,ROTATE_G2R,ROTATE_R2G}; + TRICLINIC_GENERAL,TRICLINIC_GENERAL_FLIP,QUAT_G2R}; #define LB_FACTOR 1.1 diff --git a/src/read_restart.cpp b/src/read_restart.cpp index 372d1bcfe8..ca5eb25762 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -28,6 +28,7 @@ #include "improper.h" #include "irregular.h" #include "label_map.h" +#include "math_extra.h" #include "memory.h" #include "modify.h" #include "pair.h" @@ -136,6 +137,21 @@ 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(); @@ -785,10 +801,10 @@ void ReadRestart::header() } else if (flag == TRICLINIC_GENERAL) { domain->triclinic_general = read_int(); - } else if (flag == ROTATE_G2R) { - read_double_vec(9,&domain->rotate_g2r[0][0]); - } else if (flag == ROTATE_R2G) { - read_double_vec(9,&domain->rotate_r2g[0][0]); + } 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 == SPECIAL_LJ) { read_int(); diff --git a/src/write_restart.cpp b/src/write_restart.cpp index fb5b41eb02..16a59ebc20 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -450,8 +450,8 @@ void WriteRestart::header() write_int(TRICLINIC_GENERAL,domain->triclinic_general); if (domain->triclinic_general) { - write_double_vec(ROTATE_G2R,9,&domain->rotate_g2r[0][0]); - write_double_vec(ROTATE_R2G,9,&domain->rotate_r2g[0][0]); + write_int(TRICLINIC_GENERAL_FLIP,domain->triclinic_general_flip); + write_double_vec(QUAT_G2R,4,domain->quat_g2r); } write_double_vec(SPECIAL_LJ,3,&force->special_lj[1]);