more updates

This commit is contained in:
Steve Plimpton
2023-10-12 06:49:59 -06:00
parent 94a0d85b33
commit 201f8cda9a
8 changed files with 39 additions and 225 deletions

View File

@ -53,8 +53,6 @@ AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp)
nlocal_bonus = nghost_bonus = nmax_bonus = 0; nlocal_bonus = nghost_bonus = nmax_bonus = 0;
bonus = nullptr; bonus = nullptr;
quat_hold = nullptr;
bptr = nullptr; bptr = nullptr;
if (sizeof(double) == sizeof(int)) if (sizeof(double) == sizeof(int))
@ -551,90 +549,6 @@ void AtomVecBody::data_atom_post(int ilocal)
angmom[ilocal][2] = 0.0; 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 unpack one body from Bodies section of data file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

View File

@ -60,9 +60,6 @@ class AtomVecBody : public AtomVec {
void create_atom_post(int) override; void create_atom_post(int) override;
void data_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_pre(int) override;
void pack_data_post(int) override; void pack_data_post(int) override;
@ -81,8 +78,6 @@ class AtomVecBody : public AtomVec {
double *rmass, *radius; double *rmass, *radius;
double **angmom; double **angmom;
double **quat_hold;
int nghost_bonus, nmax_bonus; int nghost_bonus, nmax_bonus;
int intdoubleratio; // sizeof(double) / sizeof(int) int intdoubleratio; // sizeof(double) / sizeof(int)
int body_flag; int body_flag;

View File

@ -49,8 +49,6 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp) : AtomVec(lmp)
nlocal_bonus = nghost_bonus = nmax_bonus = 0; nlocal_bonus = nghost_bonus = nmax_bonus = 0;
bonus = nullptr; bonus = nullptr;
quat_hold = nullptr;
// strings with peratom variables to include in each AtomVec method // strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings // strings cannot contain fields in corresponding AtomVec default strings
// order of fields in a string does not matter // 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; 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 modify values for AtomVec::pack_data() to pack
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

View File

@ -53,9 +53,6 @@ class AtomVecEllipsoid : public AtomVec {
void create_atom_post(int) override; void create_atom_post(int) override;
void data_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_pre(int) override;
void pack_data_post(int) override; void pack_data_post(int) override;
@ -73,8 +70,6 @@ class AtomVecEllipsoid : public AtomVec {
double *rmass; double *rmass;
double **angmom; double **angmom;
double **quat_hold;
int nghost_bonus, nmax_bonus; int nghost_bonus, nmax_bonus;
int ellipsoid_flag; int ellipsoid_flag;
double rmass_one; double rmass_one;

View File

@ -632,17 +632,17 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller,
if (bvec1[2] > 0.0) theta2 = -theta2; if (bvec1[2] > 0.0) theta2 = -theta2;
MathExtra::axisangle_to_quat(xaxis,theta2,quat2); 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 // quat_r2g = rotation from restricted to general
// rotate_g2r = general to restricted rotation matrix // rotate_g2r = general to restricted rotation matrix
// include flip of C vector if needed to obey right-hand rule // include flip of C vector if needed to obey right-hand rule
// rotate_r2g = restricted to general rotation matrix // rotate_r2g = restricted to general rotation matrix
// simply a transpose of rotate_g2r since orthonormal // simply a transpose of rotate_g2r since orthonormal
MathExtra::quatquat(quat2,quat1,quat_g2r); double quat_single[4];
MathExtra::qconjugate(quat_g2r,quat_r2g); 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) { if (triclinic_general_flip) {
rotate_g2r[2][0] = -rotate_g2r[2][0]; 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("Rotvec1: %g %g %g\n",rot1[0],rot1[1],rot1[2]);
printf("Theta2: %g\n",theta2); printf("Theta2: %g\n",theta2);
printf("Rotvec2: %g %g %g\n",xaxis[0],xaxis[1],xaxis[2]); 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]); printf("Quat: %g %g %g %g\n",
double angle = 2.0*acos(quat_g2r[0]); 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("Theta: %g\n",angle);
printf("Rotvec: %g %g %g\n",quat_g2r[1]/sin(0.5*angle),quat_g2r[2]/sin(0.5*angle), printf("Rotvec: %g %g %g\n",
quat_g2r[3]/sin(0.5*angle)); 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("Aprime: %g %g %g\n",aprime[0],aprime[1],aprime[2]);
printf("Bprime: %g %g %g\n",bprime[0],bprime[1],bprime[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]); printf("Cprime: %g %g %g\n",cprime[0],cprime[1],cprime[2]);

View File

@ -93,7 +93,6 @@ class Domain : protected Pointers {
// boxlo = lower left corner // boxlo = lower left corner
double avec[3], bvec[3], cvec[3]; // ABC edge vectors of general triclinic box 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_g2r[3][3]; // rotation matrix from general --> restricted tri
double rotate_r2g[3][3]; // rotation matrix from restricted --> general tri double rotate_r2g[3][3]; // rotation matrix from restricted --> general tri

View File

@ -91,9 +91,6 @@ void DumpAtom::init_style()
// setup function ptrs // 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) if (binary && domain->triclinic == 0)
header_choice = &DumpAtom::header_binary; header_choice = &DumpAtom::header_binary;
else if (binary && triclinic_general == 1) else if (binary && triclinic_general == 1)
@ -107,26 +104,36 @@ void DumpAtom::init_style()
else if (!binary && domain->triclinic == 1) else if (!binary && domain->triclinic == 1)
header_choice = &DumpAtom::header_item_triclinic; header_choice = &DumpAtom::header_item_triclinic;
if (scale_flag == 1 && image_flag == 0 && domain->triclinic == 0) if (scale_flag == 0) {
pack_choice = &DumpAtom::pack_scale_noimage; if (image_flag == 0) {
else if (scale_flag == 1 && image_flag == 1 && domain->triclinic == 0) if (triclinic_general == 1) {
pack_choice = &DumpAtom::pack_scale_image; 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; if (image_flag == 0) convert_choice = &DumpAtom::convert_noimage;
else convert_choice = &DumpAtom::convert_image; else convert_choice = &DumpAtom::convert_image;

View File

@ -1354,20 +1354,14 @@ int DumpCustom::parse_fields(int narg, char **arg)
else pack_choice[iarg] = &DumpCustom::pack_z; else pack_choice[iarg] = &DumpCustom::pack_z;
vtype[iarg] = Dump::DOUBLE; vtype[iarg] = Dump::DOUBLE;
} else if (strcmp(arg[iarg],"xs") == 0) { } 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; if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xs_triclinic;
else pack_choice[iarg] = &DumpCustom::pack_xs; else pack_choice[iarg] = &DumpCustom::pack_xs;
vtype[iarg] = Dump::DOUBLE; vtype[iarg] = Dump::DOUBLE;
} else if (strcmp(arg[iarg],"ys") == 0) { } 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; if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_ys_triclinic;
else pack_choice[iarg] = &DumpCustom::pack_ys; else pack_choice[iarg] = &DumpCustom::pack_ys;
vtype[iarg] = Dump::DOUBLE; vtype[iarg] = Dump::DOUBLE;
} else if (strcmp(arg[iarg],"zs") == 0) { } 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; if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_zs_triclinic;
else pack_choice[iarg] = &DumpCustom::pack_zs; else pack_choice[iarg] = &DumpCustom::pack_zs;
vtype[iarg] = Dump::DOUBLE; vtype[iarg] = Dump::DOUBLE;
@ -1387,20 +1381,14 @@ int DumpCustom::parse_fields(int narg, char **arg)
else pack_choice[iarg] = &DumpCustom::pack_zu; else pack_choice[iarg] = &DumpCustom::pack_zu;
vtype[iarg] = Dump::DOUBLE; vtype[iarg] = Dump::DOUBLE;
} else if (strcmp(arg[iarg],"xsu") == 0) { } 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; if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_xsu_triclinic;
else pack_choice[iarg] = &DumpCustom::pack_xsu; else pack_choice[iarg] = &DumpCustom::pack_xsu;
vtype[iarg] = Dump::DOUBLE; vtype[iarg] = Dump::DOUBLE;
} else if (strcmp(arg[iarg],"ysu") == 0) { } 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; if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_ysu_triclinic;
else pack_choice[iarg] = &DumpCustom::pack_ysu; else pack_choice[iarg] = &DumpCustom::pack_ysu;
vtype[iarg] = Dump::DOUBLE; vtype[iarg] = Dump::DOUBLE;
} else if (strcmp(arg[iarg],"zsu") == 0) { } 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; if (domain->triclinic) pack_choice[iarg] = &DumpCustom::pack_zsu_triclinic;
else pack_choice[iarg] = &DumpCustom::pack_zsu; else pack_choice[iarg] = &DumpCustom::pack_zsu;
vtype[iarg] = Dump::DOUBLE; vtype[iarg] = Dump::DOUBLE;