From e5f3fcbbf41e47decbce9deed926d4dde1bd2fd8 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Wed, 13 Sep 2023 13:29:37 -0600 Subject: [PATCH] more work on read_data and write_data --- src/DIELECTRIC/atom_vec_dielectric.cpp | 39 ++++++++++++ src/DIELECTRIC/atom_vec_dielectric.h | 5 ++ src/DIPOLE/atom_vec_dipole.cpp | 45 ++++++++++++- src/DIPOLE/atom_vec_dipole.h | 5 +- src/MACHDYN/atom_vec_smd.cpp | 41 ++++++++++++ src/MACHDYN/atom_vec_smd.h | 4 ++ src/SPIN/atom_vec_spin.cpp | 43 ++++++++++++- src/SPIN/atom_vec_spin.h | 5 +- src/atom_vec_body.cpp | 88 ++++++++++++++++++++++++++ src/atom_vec_body.h | 5 ++ src/atom_vec_ellipsoid.cpp | 87 +++++++++++++++++++++++++ src/atom_vec_ellipsoid.h | 5 ++ src/atom_vec_line.cpp | 23 +++++++ src/atom_vec_line.h | 1 + src/atom_vec_tri.cpp | 86 +++++++++++++++++++++++++ src/atom_vec_tri.h | 5 ++ src/domain.cpp | 14 ++-- src/domain.h | 1 + 18 files changed, 488 insertions(+), 14 deletions(-) diff --git a/src/DIELECTRIC/atom_vec_dielectric.cpp b/src/DIELECTRIC/atom_vec_dielectric.cpp index 735d770b04..516c08bd98 100644 --- a/src/DIELECTRIC/atom_vec_dielectric.cpp +++ b/src/DIELECTRIC/atom_vec_dielectric.cpp @@ -18,6 +18,7 @@ #include "domain.h" #include "error.h" #include "force.h" +#include "memory.h" #include "pair.h" #include "pair_hybrid.h" @@ -51,6 +52,8 @@ AtomVecDielectric::AtomVecDielectric(LAMMPS *_lmp) : AtomVec(_lmp) atom->molecule_flag = atom->q_flag = atom->mu_flag = 1; atom->dielectric_flag = 1; + mu_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 @@ -202,6 +205,42 @@ void AtomVecDielectric::read_data_general_to_restricted(int nlocal_previous, int domain->general_to_restricted_vector(mu[i]); } +/* ---------------------------------------------------------------------- + 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 dipole momemt mu +------------------------------------------------------------------------- */ + +void AtomVecDielectric::write_data_restricted_to_general() +{ + AtomVec::write_data_restricted_to_general(); + + int nlocal = atom->nlocal; + memory->create(mu_hold,nlocal,3,"atomvec:mu_hold"); + if (nlocal) memcpy(&mu_hold[0][0],&mu[0][0],3*nlocal*sizeof(double)); + for (int i = 0; i < nlocal; i++) + domain->restricted_to_general_vector(mu[i]); +} + +/* ---------------------------------------------------------------------- + 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 dipole moment mu +------------------------------------------------------------------------- */ + +void AtomVecDielectric::write_data_restore_restricted() +{ + AtomVec::write_data_restore_restricted(); + + if (!mu_hold) return; + + int nlocal = atom->nlocal; + memcpy(&mu[0][0],&mu_hold[0][0],3*nlocal*sizeof(double)); + memory->destroy(mu_hold); + mu_hold = nullptr; +} + /* ---------------------------------------------------------------------- initialize other atom quantities after AtomVec::unpack_restart() ------------------------------------------------------------------------- */ diff --git a/src/DIELECTRIC/atom_vec_dielectric.h b/src/DIELECTRIC/atom_vec_dielectric.h index 8bef111cb4..b6b7ebd676 100644 --- a/src/DIELECTRIC/atom_vec_dielectric.h +++ b/src/DIELECTRIC/atom_vec_dielectric.h @@ -36,6 +36,9 @@ class AtomVecDielectric : virtual 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 unpack_restart_init(int) override; int property_atom(const std::string &) override; void pack_property_atom(int, double *, int, int) override; @@ -49,6 +52,8 @@ class AtomVecDielectric : virtual public AtomVec { double **mu; double *area, *ed, *em, *epsilon, *curvature, *q_scaled; + + double **mu_hold; }; } // namespace LAMMPS_NS diff --git a/src/DIPOLE/atom_vec_dipole.cpp b/src/DIPOLE/atom_vec_dipole.cpp index 5323e33f17..025624c6c0 100644 --- a/src/DIPOLE/atom_vec_dipole.cpp +++ b/src/DIPOLE/atom_vec_dipole.cpp @@ -15,6 +15,7 @@ #include "atom.h" #include "domain.h" +#include "memory.h" #include @@ -28,6 +29,8 @@ AtomVecDipole::AtomVecDipole(LAMMPS *lmp) : AtomVec(lmp) mass_type = PER_TYPE; atom->q_flag = atom->mu_flag = 1; + + mu_hold = nullptr; // strings with peratom variables to include in each AtomVec method // strings cannot contain fields in corresponding AtomVec default strings @@ -73,13 +76,49 @@ void AtomVecDipole::data_atom_post(int ilocal) /* ---------------------------------------------------------------------- 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 mu + child class operates on dipole moment mu ------------------------------------------------------------------------- */ -void AtomVecDipole::data_general_to_restricted(int nlocal_previous, int nlocal) +void AtomVecDipole::read_data_general_to_restricted(int nlocal_previous, int nlocal) { - AtomVec::data_general_to_restricted(nlocal_previous, nlocal); + AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal); for (int i = nlocal_previous; i < nlocal; i++) domain->general_to_restricted_vector(mu[i]); } + +/* ---------------------------------------------------------------------- + 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 dipole momemt mu +------------------------------------------------------------------------- */ + +void AtomVecDipole::write_data_restricted_to_general() +{ + AtomVec::write_data_restricted_to_general(); + + int nlocal = atom->nlocal; + memory->create(mu_hold,nlocal,3,"atomvec:mu_hold"); + if (nlocal) memcpy(&mu_hold[0][0],&mu[0][0],3*nlocal*sizeof(double)); + for (int i = 0; i < nlocal; i++) + domain->restricted_to_general_vector(mu[i]); +} + +/* ---------------------------------------------------------------------- + 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 dipole moment mu +------------------------------------------------------------------------- */ + +void AtomVecDipole::write_data_restore_restricted() +{ + AtomVec::write_data_restore_restricted(); + + if (!mu_hold) return; + + int nlocal = atom->nlocal; + memcpy(&mu[0][0],&mu_hold[0][0],3*nlocal*sizeof(double)); + memory->destroy(mu_hold); + mu_hold = nullptr; +} diff --git a/src/DIPOLE/atom_vec_dipole.h b/src/DIPOLE/atom_vec_dipole.h index d688fd98dd..1f6d6fe2be 100644 --- a/src/DIPOLE/atom_vec_dipole.h +++ b/src/DIPOLE/atom_vec_dipole.h @@ -30,10 +30,13 @@ class AtomVecDipole : virtual public AtomVec { void grow_pointers() override; void data_atom_post(int) override; - void data_general_to_restricted(int, int); + void read_data_general_to_restricted(int, int); + void write_data_restricted_to_general(); + void write_data_restore_restricted(); protected: double **mu; + double **mu_hold; }; } // namespace LAMMPS_NS diff --git a/src/MACHDYN/atom_vec_smd.cpp b/src/MACHDYN/atom_vec_smd.cpp index 27f23c5362..5ba2c01038 100644 --- a/src/MACHDYN/atom_vec_smd.cpp +++ b/src/MACHDYN/atom_vec_smd.cpp @@ -25,6 +25,8 @@ #include "atom_vec_smd.h" #include "atom.h" +#include "domain.h" +#include "memory.h" #include @@ -57,6 +59,8 @@ AtomVecSMD::AtomVecSMD(LAMMPS *lmp) : AtomVec(lmp) atom->damage_flag = 1; atom->eff_plastic_strain_rate_flag = 1; + x0_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 @@ -184,3 +188,40 @@ void AtomVecSMD::data_atom_post(int ilocal) smd_data_9[ilocal][4] = 1.0; // yy smd_data_9[ilocal][8] = 1.0; // zz } + +/* ---------------------------------------------------------------------- + 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 original coords x0 +------------------------------------------------------------------------- */ + +void AtomVecSMD::write_data_restricted_to_general() +{ + AtomVec::write_data_restricted_to_general(); + + int nlocal = atom->nlocal; + memory->create(x0_hold,nlocal,3,"atomvec:x0_hold"); + if (nlocal) memcpy(&x0_hold[0][0],&x0[0][0],3*nlocal*sizeof(double)); + for (int i = 0; i < nlocal; i++) + domain->restricted_to_general_coords(x0[i]); + +} + +/* ---------------------------------------------------------------------- + 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 original coords x0 +------------------------------------------------------------------------- */ + +void AtomVecSMD::write_data_restore_restricted() +{ + AtomVec::write_data_restore_restricted(); + + if (!x0_hold) return; + + int nlocal = atom->nlocal; + memcpy(&x0[0][0],&x0_hold[0][0],3*nlocal*sizeof(double)); + memory->destroy(x0_hold); + x0_hold = nullptr; +} diff --git a/src/MACHDYN/atom_vec_smd.h b/src/MACHDYN/atom_vec_smd.h index 6ca7f08b4d..322136ebd3 100644 --- a/src/MACHDYN/atom_vec_smd.h +++ b/src/MACHDYN/atom_vec_smd.h @@ -43,12 +43,16 @@ class AtomVecSMD : virtual public AtomVec { void force_clear(int, size_t) override; void create_atom_post(int) override; void data_atom_post(int) override; + void write_data_restricted_to_general(); + void write_data_restore_restricted(); private: tagint *molecule; double *esph, *desph, *vfrac, *rmass, *radius, *contact_radius; double *eff_plastic_strain, *eff_plastic_strain_rate, *damage; double **x0, **smd_data_9, **smd_stress, **vest; + + double **x0_hold; }; } // namespace LAMMPS_NS diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp index 1ea57516f6..a68a59ef38 100644 --- a/src/SPIN/atom_vec_spin.cpp +++ b/src/SPIN/atom_vec_spin.cpp @@ -25,6 +25,7 @@ #include "atom.h" #include "domain.h" +#include "memory.h" #include #include @@ -41,6 +42,8 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp) atom->sp_flag = 1; + sp_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 @@ -106,10 +109,46 @@ void AtomVecSpin::data_atom_post(int ilocal) child class operates on spin vector sp ------------------------------------------------------------------------- */ -void AtomVecSpin::data_general_to_restricted(int nlocal_previous, int nlocal) +void AtomVecSpin::read_data_general_to_restricted(int nlocal_previous, int nlocal) { - AtomVec::data_general_to_restricted(nlocal_previous, nlocal); + AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal); for (int i = nlocal_previous; i < nlocal; i++) domain->general_to_restricted_vector(sp[i]); } + +/* ---------------------------------------------------------------------- + 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 spin vector sp +------------------------------------------------------------------------- */ + +void AtomVecSpin::write_data_restricted_to_general() +{ + AtomVec::write_data_restricted_to_general(); + + int nlocal = atom->nlocal; + memory->create(sp_hold,nlocal,3,"atomvec:sp_hold"); + if (nlocal) memcpy(&sp_hold[0][0],&sp[0][0],3*nlocal*sizeof(double)); + for (int i = 0; i < nlocal; i++) + domain->restricted_to_general_vector(sp[i]); +} + +/* ---------------------------------------------------------------------- + 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 spin vector sp +------------------------------------------------------------------------- */ + +void AtomVecSpin::write_data_restore_restricted() +{ + AtomVec::write_data_restore_restricted(); + + if (!sp_hold) return; + + int nlocal = atom->nlocal; + memcpy(&sp[0][0],&sp_hold[0][0],3*nlocal*sizeof(double)); + memory->destroy(sp_hold); + sp_hold = nullptr; +} diff --git a/src/SPIN/atom_vec_spin.h b/src/SPIN/atom_vec_spin.h index effbe232f4..93bbc82ab8 100644 --- a/src/SPIN/atom_vec_spin.h +++ b/src/SPIN/atom_vec_spin.h @@ -31,10 +31,13 @@ class AtomVecSpin : virtual public AtomVec { void grow_pointers() override; void force_clear(int, size_t) override; void data_atom_post(int) override; - void data_general_to_restricted(int, int); + void read_data_general_to_restricted(int, int); + void write_data_restricted_to_general(); + void write_data_restore_restricted(); protected: double **sp, **fm, **fm_long; + double **sp_hold; }; } // namespace LAMMPS_NS diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp index 538e9783df..7721e1540b 100644 --- a/src/atom_vec_body.cpp +++ b/src/atom_vec_body.cpp @@ -16,8 +16,10 @@ #include "atom.h" #include "body.h" +#include "domain.h" #include "error.h" #include "fix.h" +#include "math_extra.h" #include "memory.h" #include "modify.h" #include "my_pool_chunk.h" @@ -51,6 +53,8 @@ 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)) @@ -547,6 +551,90 @@ 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],quat_r2g[4]; + double *bquat; + + 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; + 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 e02fd3bbb0..955b4f4587 100644 --- a/src/atom_vec_body.h +++ b/src/atom_vec_body.h @@ -60,6 +60,9 @@ 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; @@ -78,6 +81,8 @@ 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 3eaa927384..ceb813cbe2 100644 --- a/src/atom_vec_ellipsoid.cpp +++ b/src/atom_vec_ellipsoid.cpp @@ -18,6 +18,7 @@ #include "atom_vec_ellipsoid.h" #include "atom.h" +#include "domain.h" #include "error.h" #include "fix.h" #include "math_const.h" @@ -48,6 +49,8 @@ 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 @@ -457,6 +460,90 @@ 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],quat_r2g[4]; + double *bquat; + + 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; + 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 6e06d773fc..3d6815fff0 100644 --- a/src/atom_vec_ellipsoid.h +++ b/src/atom_vec_ellipsoid.h @@ -53,6 +53,9 @@ 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; @@ -70,6 +73,8 @@ 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/atom_vec_line.cpp b/src/atom_vec_line.cpp index ff09bed6d0..95ccfe31bd 100644 --- a/src/atom_vec_line.cpp +++ b/src/atom_vec_line.cpp @@ -437,6 +437,29 @@ void AtomVecLine::data_atom_post(int ilocal) omega[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 theta +------------------------------------------------------------------------- */ + +void AtomVecLine::read_data_general_to_restricted(int nlocal_previous, int nlocal) +{ + AtomVec::read_data_general_to_restricted(nlocal_previous, nlocal); + + double btheta; + double theta_g2r = 2.0*acos(domain->quat_g2r[0]); + + for (int i = nlocal_previous; i < nlocal; i++) { + if (line[i] < 0) continue; + btheta = bonus[line[i]].theta; + btheta += theta_g2r; + if (btheta > MathConst::MY_PI) btheta -= MathConst::MY_2PI; + else if (btheta <= -MathConst::MY_PI) btheta += MathConst::MY_2PI; + bonus[line[i]].theta = btheta; + } +} + /* ---------------------------------------------------------------------- modify values for AtomVec::pack_data() to pack ------------------------------------------------------------------------- */ diff --git a/src/atom_vec_line.h b/src/atom_vec_line.h index 740c541916..2503f55d8d 100644 --- a/src/atom_vec_line.h +++ b/src/atom_vec_line.h @@ -53,6 +53,7 @@ class AtomVecLine : public AtomVec { void create_atom_post(int) override; void data_atom_post(int) override; + void read_data_general_to_restricted(int, int); void pack_data_pre(int) override; void pack_data_post(int) override; diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp index a46609b02c..7dfb50093b 100644 --- a/src/atom_vec_tri.cpp +++ b/src/atom_vec_tri.cpp @@ -52,6 +52,8 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp) nlocal_bonus = nghost_bonus = nmax_bonus = 0; bonus = nullptr; + double **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 @@ -657,6 +659,90 @@ void AtomVecTri::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 AtomVecTri::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 (tri[i] < 0) continue; + bquat = bonus[tri[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 AtomVecTri::write_data_restricted_to_general() +{ + AtomVec::write_data_restricted_to_general(); + + double quat[4],quat_r2g[4]; + double *bquat; + + 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; + 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 AtomVecTri::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_tri.h b/src/atom_vec_tri.h index 424bd8ea0a..91a7a20e6d 100644 --- a/src/atom_vec_tri.h +++ b/src/atom_vec_tri.h @@ -55,6 +55,9 @@ class AtomVecTri : 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; @@ -72,6 +75,8 @@ class AtomVecTri : public AtomVec { double *radius, *rmass; double **omega, **angmom; + double **quat_hold; + int nghost_bonus, nmax_bonus; int tri_flag; double rmass_one; diff --git a/src/domain.cpp b/src/domain.cpp index 810814a80f..9c6e1a936c 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -627,7 +627,7 @@ 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 = transformation via single quat = quat2 * quat1 + // 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, @@ -635,9 +635,8 @@ void Domain::define_general_triclinic(double *avec_caller, double *bvec_caller, // rotate_r2g = restricted to general transformation matrix // simply a transpose of rotate_g2r since orthonormal - double quat[4]; - MathExtra::quatquat(quat2,quat1,quat); - MathExtra::quat_to_mat(quat,rotate_g2r); + MathExtra::quatquat(quat2,quat1,quat_g2r); + MathExtra::quat_to_mat(quat_g2r,rotate_g2r); if (dot < 0.0) { rotate_g2r[2][0] = -rotate_g2r[2][0]; @@ -670,10 +669,11 @@ 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[0],quat[1],quat[2],quat[3]); - double angle = 2.0*acos(quat[0]); + 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("Theta: %g\n",angle); - printf("Rotvec: %g %g %g\n",quat[1]/sin(0.5*angle),quat[2]/sin(0.5*angle),quat[3]/sin(0.5*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("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 7b306a9426..3e8b215d29 100644 --- a/src/domain.h +++ b/src/domain.h @@ -92,6 +92,7 @@ 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]; // quaternion for general --> restricted rotation double rotate_g2r[3][3]; // rotation matrix from general --> restricted tri double rotate_r2g[3][3]; // rotation matrix from restricted --> general tri