From 201f8cda9a72f25bc6ad9d7da7afc0ed97c62ff8 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 12 Oct 2023 06:49:59 -0600 Subject: [PATCH] more updates --- src/atom_vec_body.cpp | 86 -------------------------------------- src/atom_vec_body.h | 5 --- src/atom_vec_ellipsoid.cpp | 86 -------------------------------------- src/atom_vec_ellipsoid.h | 5 --- src/domain.cpp | 18 ++++---- src/domain.h | 1 - src/dump_atom.cpp | 51 ++++++++++++---------- src/dump_custom.cpp | 12 ------ 8 files changed, 39 insertions(+), 225 deletions(-) diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp index 898b065749..fd590ad5fe 100644 --- a/src/atom_vec_body.cpp +++ b/src/atom_vec_body.cpp @@ -53,8 +53,6 @@ AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp) nlocal_bonus = nghost_bonus = nmax_bonus = 0; bonus = nullptr; - quat_hold = nullptr; - bptr = nullptr; if (sizeof(double) == sizeof(int)) @@ -551,90 +549,6 @@ void AtomVecBody::data_atom_post(int ilocal) angmom[ilocal][2] = 0.0; } -/* ---------------------------------------------------------------------- - convert read_data file info from general to restricted triclinic - parent class operates on data from Velocities section of data file - child class operates on bonus quat -------------------------------------------------------------------------- */ - -void AtomVecBody::read_data_general_to_restricted(int nlocal_previous, int nlocal) -{ - AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal); - - double quat[4]; - double *bquat; - - for (int i = nlocal_previous; i < nlocal; i++) { - if (body[i] < 0) continue; - bquat = bonus[body[i]].quat; - MathExtra::quatquat(domain->quat_g2r,bquat,quat); - bquat[0] = quat[0]; - bquat[1] = quat[1]; - bquat[2] = quat[2]; - bquat[3] = quat[3]; - MathExtra::qnormalize(bquat); - } -} - -/* ---------------------------------------------------------------------- - convert info output by write_data from restricted to general triclinic - parent class operates on x and data from Velocities section of data file - child class operates on bonus quat -------------------------------------------------------------------------- */ - -void AtomVecBody::write_data_restricted_to_general() -{ - AtomVec::write_data_restricted_to_general(); - - double quat[4]; - double *bquat; - double *quat_r2g = domain->quat_r2g; - - memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold"); - - for (int i = 0; i < nlocal_bonus; i++) { - bquat = bonus[i].quat; - quat_hold[i][0] = bquat[0]; - quat_hold[i][1] = bquat[1]; - quat_hold[i][2] = bquat[2]; - quat_hold[i][3] = bquat[3]; - - MathExtra::quatquat(quat_r2g,bquat,quat); - bquat[0] = quat[0]; - bquat[1] = quat[1]; - bquat[2] = quat[2]; - bquat[3] = quat[3]; - MathExtra::qnormalize(bquat); - } -} - -/* ---------------------------------------------------------------------- - restore info output by write_data to restricted triclinic - original data is in "hold" arrays - parent class operates on x and data from Velocities section of data file - child class operates on bonus quat -------------------------------------------------------------------------- */ - -void AtomVecBody::write_data_restore_restricted() -{ - AtomVec::write_data_restore_restricted(); - - if (!quat_hold) return; - - double *bquat; - - for (int i = 0; i < nlocal_bonus; i++) { - bquat = bonus[i].quat; - bquat[0] = quat_hold[i][0]; - bquat[1] = quat_hold[i][1]; - bquat[2] = quat_hold[i][2]; - bquat[3] = quat_hold[i][3]; - } - - memory->destroy(quat_hold); - quat_hold = nullptr; -} - /* ---------------------------------------------------------------------- unpack one body from Bodies section of data file ------------------------------------------------------------------------- */ diff --git a/src/atom_vec_body.h b/src/atom_vec_body.h index 955b4f4587..e02fd3bbb0 100644 --- a/src/atom_vec_body.h +++ b/src/atom_vec_body.h @@ -60,9 +60,6 @@ class AtomVecBody : public AtomVec { void create_atom_post(int) override; void data_atom_post(int) override; - void read_data_general_to_restricted(int, int); - void write_data_restricted_to_general(); - void write_data_restore_restricted(); void pack_data_pre(int) override; void pack_data_post(int) override; @@ -81,8 +78,6 @@ class AtomVecBody : public AtomVec { double *rmass, *radius; double **angmom; - double **quat_hold; - int nghost_bonus, nmax_bonus; int intdoubleratio; // sizeof(double) / sizeof(int) int body_flag; diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp index d51f7af78f..3cc8f6362d 100644 --- a/src/atom_vec_ellipsoid.cpp +++ b/src/atom_vec_ellipsoid.cpp @@ -49,8 +49,6 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp) : AtomVec(lmp) nlocal_bonus = nghost_bonus = nmax_bonus = 0; bonus = nullptr; - quat_hold = nullptr; - // strings with peratom variables to include in each AtomVec method // strings cannot contain fields in corresponding AtomVec default strings // order of fields in a string does not matter @@ -460,90 +458,6 @@ void AtomVecEllipsoid::data_atom_post(int ilocal) angmom[ilocal][2] = 0.0; } -/* ---------------------------------------------------------------------- - convert read_data file info from general to restricted triclinic - parent class operates on data from Velocities section of data file - child class operates on bonus quat -------------------------------------------------------------------------- */ - -void AtomVecEllipsoid::read_data_general_to_restricted(int nlocal_previous, int nlocal) -{ - AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal); - - double quat[4]; - double *bquat; - - for (int i = nlocal_previous; i < nlocal; i++) { - if (ellipsoid[i] < 0) continue; - bquat = bonus[ellipsoid[i]].quat; - MathExtra::quatquat(domain->quat_g2r,bquat,quat); - bquat[0] = quat[0]; - bquat[1] = quat[1]; - bquat[2] = quat[2]; - bquat[3] = quat[3]; - MathExtra::qnormalize(bquat); - } -} - -/* ---------------------------------------------------------------------- - convert info output by write_data from restricted to general triclinic - parent class operates on x and data from Velocities section of data file - child class operates on bonus quat -------------------------------------------------------------------------- */ - -void AtomVecEllipsoid::write_data_restricted_to_general() -{ - AtomVec::write_data_restricted_to_general(); - - double quat[4]; - double *bquat; - double *quat_r2g = domain->quat_r2g; - - memory->create(quat_hold,nlocal_bonus,4,"atomvec:quat_hold"); - - for (int i = 0; i < nlocal_bonus; i++) { - bquat = bonus[i].quat; - quat_hold[i][0] = bquat[0]; - quat_hold[i][1] = bquat[1]; - quat_hold[i][2] = bquat[2]; - quat_hold[i][3] = bquat[3]; - - MathExtra::quatquat(quat_r2g,bquat,quat); - bquat[0] = quat[0]; - bquat[1] = quat[1]; - bquat[2] = quat[2]; - bquat[3] = quat[3]; - MathExtra::qnormalize(bquat); - } -} - -/* ---------------------------------------------------------------------- - restore info output by write_data to restricted triclinic - original data is in "hold" arrays - parent class operates on x and data from Velocities section of data file - child class operates on bonus quat -------------------------------------------------------------------------- */ - -void AtomVecEllipsoid::write_data_restore_restricted() -{ - AtomVec::write_data_restore_restricted(); - - if (!quat_hold) return; - - double *bquat; - - for (int i = 0; i < nlocal_bonus; i++) { - bquat = bonus[i].quat; - bquat[0] = quat_hold[i][0]; - bquat[1] = quat_hold[i][1]; - bquat[2] = quat_hold[i][2]; - bquat[3] = quat_hold[i][3]; - } - - memory->destroy(quat_hold); - quat_hold = nullptr; -} - /* ---------------------------------------------------------------------- modify values for AtomVec::pack_data() to pack ------------------------------------------------------------------------- */ diff --git a/src/atom_vec_ellipsoid.h b/src/atom_vec_ellipsoid.h index 3d6815fff0..6e06d773fc 100644 --- a/src/atom_vec_ellipsoid.h +++ b/src/atom_vec_ellipsoid.h @@ -53,9 +53,6 @@ class AtomVecEllipsoid : public AtomVec { void create_atom_post(int) override; void data_atom_post(int) override; - void read_data_general_to_restricted(int, int); - void write_data_restricted_to_general(); - void write_data_restore_restricted(); void pack_data_pre(int) override; void pack_data_post(int) override; @@ -73,8 +70,6 @@ class AtomVecEllipsoid : public AtomVec { double *rmass; double **angmom; - double **quat_hold; - int nghost_bonus, nmax_bonus; int ellipsoid_flag; double rmass_one; diff --git a/src/domain.cpp b/src/domain.cpp index 68eafab356..379b626d80 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -632,17 +632,17 @@ 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 = rotation via single quat = quat2 * quat1 + // quat_single = 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); + double quat_single[4]; + MathExtra::quatquat(quat2,quat1,quat_single); - MathExtra::quat_to_mat(quat_g2r,rotate_g2r); + MathExtra::quat_to_mat(quat_single,rotate_g2r); if (triclinic_general_flip) { rotate_g2r[2][0] = -rotate_g2r[2][0]; @@ -675,11 +675,13 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller, printf("Rotvec1: %g %g %g\n",rot1[0],rot1[1],rot1[2]); printf("Theta2: %g\n",theta2); printf("Rotvec2: %g %g %g\n",xaxis[0],xaxis[1],xaxis[2]); - printf("Quat: %g %g %g %g\n",quat_g2r[0],quat_g2r[1],quat_g2r[2],quat_g2r[3]); - double angle = 2.0*acos(quat_g2r[0]); + printf("Quat: %g %g %g %g\n", + quat_single[0],quat_single[1],quat_single[2],quat_single[3]); + double angle = 2.0*acos(quat_single[0]); printf("Theta: %g\n",angle); - printf("Rotvec: %g %g %g\n",quat_g2r[1]/sin(0.5*angle),quat_g2r[2]/sin(0.5*angle), - quat_g2r[3]/sin(0.5*angle)); + printf("Rotvec: %g %g %g\n", + quat_single[1]/sin(0.5*angle),quat_single[2]/sin(0.5*angle), + quat_single[3]/sin(0.5*angle)); printf("Aprime: %g %g %g\n",aprime[0],aprime[1],aprime[2]); printf("Bprime: %g %g %g\n",bprime[0],bprime[1],bprime[2]); printf("Cprime: %g %g %g\n",cprime[0],cprime[1],cprime[2]); diff --git a/src/domain.h b/src/domain.h index 61a88724e7..34714cf6e3 100644 --- a/src/domain.h +++ b/src/domain.h @@ -93,7 +93,6 @@ class Domain : protected Pointers { // boxlo = lower left corner double avec[3], bvec[3], cvec[3]; // ABC edge vectors of general triclinic box - 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 151446bdfd..38c8431269 100644 --- a/src/dump_atom.cpp +++ b/src/dump_atom.cpp @@ -91,9 +91,6 @@ void DumpAtom::init_style() // setup function ptrs - if (scale_flag && triclinic_general) - error->all(FLERR,"Dump atom cannot use scale and triclinic/general settings"); - if (binary && domain->triclinic == 0) header_choice = &DumpAtom::header_binary; else if (binary && triclinic_general == 1) @@ -107,26 +104,36 @@ void DumpAtom::init_style() else if (!binary && domain->triclinic == 1) header_choice = &DumpAtom::header_item_triclinic; - if (scale_flag == 1 && image_flag == 0 && domain->triclinic == 0) - pack_choice = &DumpAtom::pack_scale_noimage; - else if (scale_flag == 1 && image_flag == 1 && domain->triclinic == 0) - pack_choice = &DumpAtom::pack_scale_image; + if (scale_flag == 0) { + if (image_flag == 0) { + if (triclinic_general == 1) { + pack_choice = &DumpAtom::pack_noscale_noimage_triclinic_general; + } else { + pack_choice = &DumpAtom::pack_noscale_noimage; + } + } else if (image_flag == 1) { + if (triclinic_general == 1) { + pack_choice = &DumpAtom::pack_noscale_image_triclinic_general; + } else { + pack_choice = &DumpAtom::pack_noscale_image; + } + } + } else if (scale_flag == 1) { + if (image_flag == 0) { + if (domain->triclinic == 0) { + pack_choice = &DumpAtom::pack_scale_noimage; + } else { + pack_choice = &DumpAtom::pack_scale_noimage_triclinic; + } + } else if (image_flag == 1) { + if (domain->triclinic == 0) { + pack_choice = &DumpAtom::pack_scale_image; + } else { + pack_choice = &DumpAtom::pack_scale_image_triclinic; + } + } + } - else if (scale_flag == 0 && image_flag == 0 && triclinic_general == 1) - pack_choice = &DumpAtom::pack_noscale_noimage_triclinic_general; - else if (scale_flag == 0 && image_flag == 1 && triclinic_general == 1) - pack_choice = &DumpAtom::pack_noscale_image_triclinic_general; - - else if (scale_flag == 1 && image_flag == 0 && domain->triclinic == 1) - pack_choice = &DumpAtom::pack_scale_noimage_triclinic; - else if (scale_flag == 1 && image_flag == 1 && domain->triclinic == 1) - pack_choice = &DumpAtom::pack_scale_image_triclinic; - - else if (scale_flag == 0 && image_flag == 0) - pack_choice = &DumpAtom::pack_noscale_noimage; - else if (scale_flag == 0 && image_flag == 1) - pack_choice = &DumpAtom::pack_noscale_image; - if (image_flag == 0) convert_choice = &DumpAtom::convert_noimage; else convert_choice = &DumpAtom::convert_image; diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 2e478fd5e9..9554247aff 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -1354,20 +1354,14 @@ int DumpCustom::parse_fields(int narg, char **arg) else pack_choice[iarg] = &DumpCustom::pack_z; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"xs") == 0) { - if (triclinic_general) - error->all(FLERR,"Dump custom xs property not supported for general triclinic"); if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xs_triclinic; else pack_choice[iarg] = &DumpCustom::pack_xs; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"ys") == 0) { - if (triclinic_general) - error->all(FLERR,"Dump custom ys property not supported for general triclinic"); if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_ys_triclinic; else pack_choice[iarg] = &DumpCustom::pack_ys; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"zs") == 0) { - if (triclinic_general) - error->all(FLERR,"Dump custom zs property not supported for general triclinic"); if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_zs_triclinic; else pack_choice[iarg] = &DumpCustom::pack_zs; vtype[iarg] = Dump::DOUBLE; @@ -1387,20 +1381,14 @@ int DumpCustom::parse_fields(int narg, char **arg) else pack_choice[iarg] = &DumpCustom::pack_zu; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"xsu") == 0) { - if (triclinic_general) - error->all(FLERR,"Dump custom xsu property not supported for general triclinic"); if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xsu_triclinic; else pack_choice[iarg] = &DumpCustom::pack_xsu; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"ysu") == 0) { - if (triclinic_general) - error->all(FLERR,"Dump custom ysu property not supported for general triclinic"); if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_ysu_triclinic; else pack_choice[iarg] = &DumpCustom::pack_ysu; vtype[iarg] = Dump::DOUBLE; } else if (strcmp(arg[iarg],"zsu") == 0) { - if (triclinic_general) - error->all(FLERR,"Dump custom zsu property not supported for general triclinic"); if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_zsu_triclinic; else pack_choice[iarg] = &DumpCustom::pack_zsu; vtype[iarg] = Dump::DOUBLE;