initial support for write_data

This commit is contained in:
Steve Plimpton
2023-09-13 10:03:05 -06:00
parent dec245c67b
commit c7e794146f
8 changed files with 124 additions and 32 deletions

View File

@ -194,9 +194,9 @@ void AtomVecDielectric::data_atom_post(int ilocal)
child class operates on dipole moment mu child class operates on dipole moment mu
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void AtomVecDielectric::data_general_to_restricted(int nlocal_previous, int nlocal) void AtomVecDielectric::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++) for (int i = nlocal_previous; i < nlocal; i++)
domain->general_to_restricted_vector(mu[i]); domain->general_to_restricted_vector(mu[i]);

View File

@ -35,7 +35,7 @@ class AtomVecDielectric : virtual public AtomVec {
void grow_pointers() override; void grow_pointers() override;
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 data_general_to_restricted(int, int); void read_data_general_to_restricted(int, int);
void unpack_restart_init(int) override; void unpack_restart_init(int) override;
int property_atom(const std::string &) override; int property_atom(const std::string &) override;
void pack_property_atom(int, double *, int, int) override; void pack_property_atom(int, double *, int, int) override;

View File

@ -160,7 +160,7 @@ void AtomVecSMD::data_atom_post(int ilocal)
// x and x0 are in Atoms section of data file // x and x0 are in Atoms section of data file
// reset x0 b/c x may have been modified in Atom::data_atoms() // reset x0 b/c x may have been modified in Atom::data_atoms()
// for PBC, shift, etc // for PBC, shift, etc
// this also means no need for data_general_to_restricted() method // this also means no need for read_data_general_to_restricted() method
// to rotate x0 for general triclinic // to rotate x0 for general triclinic
x0[ilocal][0] = x[ilocal][0]; x0[ilocal][0] = x[ilocal][0];

View File

@ -1191,6 +1191,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
} }
// convert atom coords from general triclinic to restricted triclinic // convert atom coords from general triclinic to restricted triclinic
// so can decide which proc owns the atom
if (triclinic_general) domain->general_to_restricted_coords(xdata); if (triclinic_general) domain->general_to_restricted_coords(xdata);
@ -1216,6 +1217,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
if (coord[0] >= sublo[0] && coord[0] < subhi[0] && if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] && coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) { coord[2] >= sublo[2] && coord[2] < subhi[2]) {
avec->data_atom(xdata,imagedata,values,typestr); avec->data_atom(xdata,imagedata,values,typestr);
typestr = utils::utf8_subst(typestr); typestr = utils::utf8_subst(typestr);
if (id_offset) tag[nlocal-1] += id_offset; if (id_offset) tag[nlocal-1] += id_offset;

View File

@ -68,6 +68,9 @@ AtomVec::AtomVec(LAMMPS *lmp) : Pointers(lmp)
image = nullptr; image = nullptr;
x = v = f = nullptr; x = v = f = nullptr;
x_hold = nullptr;
v_hold = omega_hold = angmom_hold = nullptr;
threads = nullptr; threads = nullptr;
} }
@ -2223,12 +2226,12 @@ void AtomVec::write_improper(FILE *fp, int n, tagint **buf, int index)
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
convert read_data file info from general to restricted triclinic convert info input by read_data from general to restricted triclinic
parent class only operates on data from Velocities section of data file parent class only operates on data from Velocities section of data file
child classes operate on all other data: Atoms, Ellipsoids, Lines, Triangles, etc child classes operate on all other data: Atoms, Ellipsoids, Lines, Triangles, etc
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void AtomVec::data_general_to_restricted(int nlocal_previous, int nlocal) void AtomVec::read_data_general_to_restricted(int nlocal_previous, int nlocal)
{ {
int datatype, cols; int datatype, cols;
void *pdata; void *pdata;
@ -2239,7 +2242,7 @@ void AtomVec::data_general_to_restricted(int nlocal_previous, int nlocal)
cols = mdata_vel.cols[n]; cols = mdata_vel.cols[n];
// operate on v, omega, angmom // operate on v, omega, angmom
// no other read_data atom fields are Nx3 double arrays // no other read_data Velocities fields are Nx3 double arrays
if (datatype == Atom::DOUBLE) { if (datatype == Atom::DOUBLE) {
if (cols == 3) { if (cols == 3) {
@ -2251,6 +2254,99 @@ void AtomVec::data_general_to_restricted(int nlocal_previous, int nlocal)
} }
} }
/* ----------------------------------------------------------------------
convert info output by write_data from restricted to general triclinic
parent class only operates on x and data from Velocities section of data file
child classes operate on all other data: Atoms, Ellipsoids, Lines, Triangles, etc
------------------------------------------------------------------------- */
void AtomVec::write_data_restricted_to_general()
{
int datatype, cols;
void *pdata;
int nlocal = atom->nlocal;
memory->create(x_hold,nlocal,3,"atomvec:x_hold");
if (nlocal) memcpy(&x_hold[0][0],&x[0][0],3*nlocal*sizeof(double));
for (int i = 0; i < nlocal; i++)
domain->restricted_to_general_coords(x[i]);
double **omega = atom->omega;
double **angmom = atom->angmom;
for (int n = 1; n < ndata_vel; n++) {
pdata = mdata_vel.pdata[n];
datatype = mdata_vel.datatype[n];
cols = mdata_vel.cols[n];
// operate on v, omega, angmom
// no other write_data Velocities fields are Nx3 double arrays
if (datatype == Atom::DOUBLE) {
if (cols == 3) {
double **array = *((double ***) pdata);
if (array == v) {
memory->create(v_hold,nlocal,3,"atomvec:v_hold");
if (nlocal) memcpy(&v_hold[0][0],&v[0][0],3*nlocal*sizeof(double));
for (int i = 0; i < nlocal; i++)
domain->restricted_to_general_vector(v[i]);
} else if (array == omega) {
memory->create(omega_hold,nlocal,3,"atomvec:omega_hold");
if (nlocal) memcpy(&omega_hold[0][0],&omega[0][0],3*nlocal*sizeof(double));
for (int i = 0; i < nlocal; i++)
domain->restricted_to_general_vector(omega[i]);
} else if (array == angmom) {
memory->create(angmom_hold,nlocal,3,"atomvec:angmom_hold");
if (nlocal) memcpy(&angmom_hold[0][0],&angmom[0][0],3*nlocal*sizeof(double));
for (int i = 0; i < nlocal; i++)
domain->restricted_to_general_vector(angmom[i]);
}
}
}
}
}
/* ----------------------------------------------------------------------
restore info output by write_data to restricted triclinic
original data is in "hold" arrays
parent class only operates on x and data from Velocities section of data file
child classes operate on all other data: Atoms, Ellipsoids, Lines, Triangles, etc
------------------------------------------------------------------------- */
void AtomVec::write_data_restore_restricted()
{
int nlocal = atom->nlocal;
if (x_hold) {
memcpy(&x[0][0],&x_hold[0][0],3*nlocal*sizeof(double));
memory->destroy(x_hold);
x_hold = nullptr;
}
// operate on v, omega, angmom
// no other write_data Velocities fields are Nx3 double arrays
if (v_hold) {
memcpy(&v[0][0],&v_hold[0][0],3*nlocal*sizeof(double));
memory->destroy(v_hold);
v_hold = nullptr;
}
if (omega_hold) {
memcpy(&atom->omega[0][0],&omega_hold[0][0],3*nlocal*sizeof(double));
memory->destroy(omega_hold);
omega_hold = nullptr;
}
if (angmom_hold) {
memcpy(&atom->angmom[0][0],&angmom_hold[0][0],3*nlocal*sizeof(double));
memory->destroy(angmom_hold);
angmom_hold = nullptr;
}
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
return # of bytes of allocated memory return # of bytes of allocated memory
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

View File

@ -152,7 +152,9 @@ class AtomVec : protected Pointers {
virtual int pack_data_bonus(double *, int) { return 0; } virtual int pack_data_bonus(double *, int) { return 0; }
virtual void write_data_bonus(FILE *, int, double *, int) {} virtual void write_data_bonus(FILE *, int, double *, int) {}
virtual void data_general_to_restricted(int, int); virtual void read_data_general_to_restricted(int, int);
virtual void write_data_restricted_to_general();
virtual void write_data_restore_restricted();
virtual int property_atom(const std::string &) { return -1; } virtual int property_atom(const std::string &) { return -1; }
virtual void pack_property_atom(int, double *, int, int) {} virtual void pack_property_atom(int, double *, int, int) {}
@ -171,6 +173,11 @@ class AtomVec : protected Pointers {
imageint *image; imageint *image;
double **x, **v, **f; double **x, **v, **f;
// copies of original unrotated fields for write_data for general triclinic
double **x_hold;
double **v_hold, **omega_hold, **angmom_hold;
// standard list of peratom fields always operated on by different methods // standard list of peratom fields always operated on by different methods
// common to all styles, so not listed in field strings // common to all styles, so not listed in field strings

View File

@ -1066,7 +1066,7 @@ void ReadData::command(int narg, char **arg)
// on any quantities read from data file which require it // on any quantities read from data file which require it
if (triclinic_general) if (triclinic_general)
atom->avec->data_general_to_restricted(nlocal_previous, atom->nlocal); atom->avec->read_data_general_to_restricted(nlocal_previous, atom->nlocal);
// init per-atom fix/compute/variable values for created atoms // init per-atom fix/compute/variable values for created atoms

View File

@ -97,7 +97,7 @@ void WriteData::command(int narg, char **arg)
} else if (strcmp(arg[iarg],"nofix") == 0) { } else if (strcmp(arg[iarg],"nofix") == 0) {
fixflag = 0; fixflag = 0;
iarg++; iarg++;
} else if (strcmp(arg[iarg],"triclinic") == 0) { } else if (strcmp(arg[iarg],"triclinic/general") == 0) {
triclinic_general = 1; triclinic_general = 1;
iarg++; iarg++;
} else if (strcmp(arg[iarg],"nolabelmap") == 0) { } else if (strcmp(arg[iarg],"nolabelmap") == 0) {
@ -213,32 +213,14 @@ void WriteData::write(const std::string &file)
if (coeffflag) force_fields(); if (coeffflag) force_fields();
} }
// per atom info in Atoms and Velocities sections
// if general triclinic: // if general triclinic:
// save restricted triclinic atom coords // reset internal per-atom data that needs rotation
// transform atom coords from restricted to general
// restore saved atom coords after output atom->avec->write_data_restricted_to_general();
double **xstore = nullptr; // per atom info in Atoms and Velocities sections
if (triclinic_general) {
double **x = atom->x;
int nlocal = atom->nlocal;
memory->create(xstore,nlocal,3,"write_data:xstore");
if (nlocal) memcpy(&xstore[0][0],&x[0][0],3*nlocal*sizeof(double));
for (int i = 0; i < nlocal; i++)
domain->restricted_to_general_coords(x[i]);
}
if (natoms) atoms(); if (natoms) atoms();
if (triclinic_general) {
double **x = atom->x;
int nlocal = atom->nlocal;
if (nlocal) memcpy(&x[0][0],&xstore[0][0],3*nlocal*sizeof(double));
memory->destroy(xstore);
}
if (natoms) velocities(); if (natoms) velocities();
// molecular topology info if defined // molecular topology info if defined
@ -265,6 +247,11 @@ void WriteData::write(const std::string &file)
if (ifix->wd_section) if (ifix->wd_section)
for (int m = 0; m < ifix->wd_section; m++) fix(ifix,m); for (int m = 0; m < ifix->wd_section; m++) fix(ifix,m);
// if general triclinic:
// restore internal per-atom data that was rotated
atom->avec->write_data_restore_restricted();
// close data file // close data file
if (me == 0) fclose(fp); if (me == 0) fclose(fp);