Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/md/lammps/atom_vec_lm.hh
43a54ac855b9master

Last Change
Raw File

atom_vec_lm.hh
View Options

/**
* @file atom_vec_lm.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @brief Atomic vector allowing to add libmultiscale bonus to atom packs
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_ATOM_VEC_LM_HH__
#define __LIBMULTISCALE_ATOM_VEC_LM_HH__
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "reference_manager_lammps.hh"
/* -------------------------------------------------------------------------- */
#include "atom_vec.h"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename AtomVecLAMMPS, UInt Dim>
class AtomVecLM : public AtomVecLAMMPS {
public:
struct Bonus {
double shape[3];
double quat[4];
int ilocal;
};
struct Bonus *bonus;
AtomVecLM(const AtomVecLAMMPS &);
~AtomVecLM();
void grow_pointers();
void copy_bonus(int, int, int);
void clear_bonus();
int pack_comm_bonus(int, int *, double *);
void unpack_comm_bonus(int, int, double *);
int pack_border_bonus(int, int *, double *);
int unpack_border_bonus(int, int, double *);
int pack_exchange_bonus(int, double *);
int unpack_exchange_bonus(int, double *);
int size_restart_bonus();
int pack_restart_bonus(int, LAMMPS_NS::PackBuffer);
int unpack_restart_bonus(int, LAMMPS_NS::PackBuffer);
void data_atom_bonus(int, char **);
LAMMPS_NS::bigint memory_usage_bonus();
void create_atom_post(int);
void data_atom_post(int);
void pack_data_pre(int);
void pack_data_post(int);
int pack_data_bonus(double *, int);
void write_data_bonus(FILE *, int, double *, int);
// unique to AtomVecEllipsoid
void set_shape(int, double, double, double);
int nlocal_bonus;
private:
int *ellipsoid;
double *rmass;
double **angmom;
int nghost_bonus, nmax_bonus;
int ellipsoid_flag;
double rmass_one;
void grow_bonus();
void copy_bonus_all(int, int);
//! shared pointer to the reference manager
std::shared_ptr<ReferenceManagerLammps<Dim>> ref_manager;
};
/* ---------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
AtomVecLM<AtomVecLAMMPS, Dim>::AtomVecLM(const AtomVecLAMMPS &avec)
: AtomVecLAMMPS(avec) {}
/* ---------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
AtomVecLM<AtomVecLAMMPS, Dim>::~AtomVecLM() {
this->memory->sfree(bonus);
}
/* ----------------------------------------------------------------------
set local copies of all grow ptrs used by this class, except defaults
needed in replicate when 2 atom classes exist and it calls pack_restart()
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::grow_pointers() {
ellipsoid = this->atom->ellipsoid;
rmass = this->atom->rmass;
angmom = this->atom->angmom;
}
/* ----------------------------------------------------------------------
grow bonus data structure
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::grow_bonus() {
nmax_bonus = this->grow_nmax_bonus(nmax_bonus);
if (nmax_bonus < 0)
this->error->one(FLERR, "Per-processor system is too big");
bonus = (Bonus *)this->memory->srealloc(bonus, nmax_bonus * sizeof(Bonus),
"atom:bonus");
}
/* ----------------------------------------------------------------------
copy atom I bonus info to atom J
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::copy_bonus(int i, int j, int delflag) {
// if deleting atom J via delflag and J has bonus data, then delete it
if (delflag && ellipsoid[j] >= 0) {
copy_bonus_all(nlocal_bonus - 1, ellipsoid[j]);
nlocal_bonus--;
}
// if atom I has bonus data, reset I's bonus.ilocal to loc J
// do NOT do this if self-copy (I=J) since I's bonus data is already deleted
if (ellipsoid[i] >= 0 && i != j)
bonus[ellipsoid[i]].ilocal = j;
ellipsoid[j] = ellipsoid[i];
}
/* ----------------------------------------------------------------------
copy bonus data from I to J, effectively deleting the J entry
also reset ellipsoid that points to I to now point to J
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::copy_bonus_all(int i, int j) {
ellipsoid[bonus[i].ilocal] = j;
memcpy(&bonus[j], &bonus[i], sizeof(Bonus));
}
/* ----------------------------------------------------------------------
clear ghost info in bonus data
called before ghosts are recommunicated in comm and irregular
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::clear_bonus() {
nghost_bonus = 0;
if (this->atom->nextra_grow)
for (int iextra = 0; iextra < this->atom->nextra_grow; iextra++)
this->modify->fix[this->atom->extra_grow[iextra]]->clear_bonus();
}
/* ---------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
int AtomVecLM<AtomVecLAMMPS, Dim>::pack_comm_bonus(int n, int *list,
double *buf) {
int i, j, m;
double *quat;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
if (ellipsoid[j] >= 0) {
quat = bonus[ellipsoid[j]].quat;
buf[m++] = quat[0];
buf[m++] = quat[1];
buf[m++] = quat[2];
buf[m++] = quat[3];
}
}
return m;
}
/* ---------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::unpack_comm_bonus(int n, int first,
double *buf) {
int i, m, last;
double *quat;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (ellipsoid[i] >= 0) {
quat = bonus[ellipsoid[i]].quat;
quat[0] = buf[m++];
quat[1] = buf[m++];
quat[2] = buf[m++];
quat[3] = buf[m++];
}
}
}
/* ---------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
int AtomVecLM<AtomVecLAMMPS, Dim>::pack_border_bonus(int n, int *list,
double *buf) {
// int i, j, m;
// double *shape, *quat;
throw;
// m = 0;
// for (i = 0; i < n; i++) {
// j = list[i];
// if (ellipsoid[j] < 0)
// buf[m++] = this->ubuf(0).d;
// else {
// buf[m++] = this->ubuf(1).d;
// shape = bonus[ellipsoid[j]].shape;
// quat = bonus[ellipsoid[j]].quat;
// buf[m++] = shape[0];
// buf[m++] = shape[1];
// buf[m++] = shape[2];
// buf[m++] = quat[0];
// buf[m++] = quat[1];
// buf[m++] = quat[2];
// buf[m++] = quat[3];
// }
// }
// return m;
}
/* ---------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
int AtomVecLM<AtomVecLAMMPS, Dim>::unpack_border_bonus(int n, int first,
double *buf) {
throw;
// int i, j, m, last;
// double *shape, *quat;
// m = 0;
// last = first + n;
// for (i = first; i < last; i++) {
// ellipsoid[i] = (int)this->ubuf(buf[m++]).i;
// if (ellipsoid[i] == 0)
// ellipsoid[i] = -1;
// else {
// j = nlocal_bonus + nghost_bonus;
// if (j == nmax_bonus)
// grow_bonus();
// shape = bonus[j].shape;
// quat = bonus[j].quat;
// shape[0] = buf[m++];
// shape[1] = buf[m++];
// shape[2] = buf[m++];
// quat[0] = buf[m++];
// quat[1] = buf[m++];
// quat[2] = buf[m++];
// quat[3] = buf[m++];
// bonus[j].ilocal = i;
// ellipsoid[i] = j;
// nghost_bonus++;
// }
// }
// return m;
}
/* ----------------------------------------------------------------------
pack data for atom I for sending to another proc
xyz must be 1st 3 values, so comm::exchange() can test on them
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
int AtomVecLM<AtomVecLAMMPS, Dim>::pack_exchange_bonus(int i, double *buf) {
int m = AtomVecLAMMPS::pack_exchange_bonus(i, buf);
m += ref_manager->pack_subset_data(i, buf);
return m;
}
/* ---------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
int AtomVecLM<AtomVecLAMMPS, Dim>::unpack_exchange_bonus(int ilocal,
double *buf) {
throw;
// int m = 0;
// ellipsoid[ilocal] = (int)this->ubuf(buf[m++]).i;
// if (ellipsoid[ilocal] == 0)
// ellipsoid[ilocal] = -1;
// else {
// if (nlocal_bonus == nmax_bonus)
// grow_bonus();
// double *shape = bonus[nlocal_bonus].shape;
// double *quat = bonus[nlocal_bonus].quat;
// shape[0] = buf[m++];
// shape[1] = buf[m++];
// shape[2] = buf[m++];
// quat[0] = buf[m++];
// quat[1] = buf[m++];
// quat[2] = buf[m++];
// quat[3] = buf[m++];
// bonus[nlocal_bonus].ilocal = ilocal;
// ellipsoid[ilocal] = nlocal_bonus++;
// }
// return m;
}
/* ----------------------------------------------------------------------
size of restart data for all atoms owned by this proc
include extra data stored by fixes
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
int AtomVecLM<AtomVecLAMMPS, Dim>::size_restart_bonus() {
int i;
int n = 0;
int nlocal = this->atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (ellipsoid[i] >= 0)
n += this->size_restart_bonus_one;
else
n++;
}
return n;
}
/* ----------------------------------------------------------------------
pack atom I's data for restart file including bonus data
xyz must be 1st 3 values, so that read_restart can test on them
molecular types may be negative, but write as positive
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
int AtomVecLM<AtomVecLAMMPS, Dim>::pack_restart_bonus(
int i, LAMMPS_NS::PackBuffer buf) {
throw;
// int m = 0;
// if (ellipsoid[i] < 0)
// buf[m++] = this->ubuf(0).d;
// else {
// buf[m++] = this->ubuf(1).d;
// int j = ellipsoid[i];
// buf[m++] = bonus[j].shape[0];
// buf[m++] = bonus[j].shape[1];
// buf[m++] = bonus[j].shape[2];
// buf[m++] = bonus[j].quat[0];
// buf[m++] = bonus[j].quat[1];
// buf[m++] = bonus[j].quat[2];
// buf[m++] = bonus[j].quat[3];
// }
// return m;
}
/* ----------------------------------------------------------------------
unpack data for one atom from restart file including bonus data
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
int AtomVecLM<AtomVecLAMMPS, Dim>::unpack_restart_bonus(
int ilocal, LAMMPS_NS::PackBuffer buf) {
throw;
// int m = 0;
// ellipsoid[ilocal] = (int)this->ubuf(buf[m++]).i;
// if (ellipsoid[ilocal] == 0)
// ellipsoid[ilocal] = -1;
// else {
// if (nlocal_bonus == nmax_bonus)
// grow_bonus();
// double *shape = bonus[nlocal_bonus].shape;
// double *quat = bonus[nlocal_bonus].quat;
// shape[0] = buf[m++];
// shape[1] = buf[m++];
// shape[2] = buf[m++];
// quat[0] = buf[m++];
// quat[1] = buf[m++];
// quat[2] = buf[m++];
// quat[3] = buf[m++];
// bonus[nlocal_bonus].ilocal = ilocal;
// ellipsoid[ilocal] = nlocal_bonus++;
// }
// return m;
}
/* ----------------------------------------------------------------------
unpack one line from Ellipsoids section of data file
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::data_atom_bonus(int m, char **) {
if (ellipsoid[m])
this->error->one(FLERR,
"Assigning ellipsoid parameters to non-ellipsoid atom");
if (nlocal_bonus == nmax_bonus)
grow_bonus();
// double *shape = bonus[nlocal_bonus].shape;
// shape[0] = 0.5 * LAMMPS_NS::utils::numeric(FLERR, values[0], true,
// this->lmp); shape[1] = 0.5 * LAMMPS_NS::utils::numeric(FLERR, values[1],
// true, this->lmp); shape[2] = 0.5 * LAMMPS_NS::utils::numeric(FLERR,
// values[2], true, this->lmp); if (shape[0] <= 0.0 || shape[1] <= 0.0 ||
// shape[2] <= 0.0)
// this->error->one(FLERR, "Invalid shape in Ellipsoids section of data
// file");
// double *quat = bonus[nlocal_bonus].quat;
// quat[0] = LAMMPS_NS::utils::numeric(FLERR, values[3], true, lmp);
// quat[1] = LAMMPS_NS::utils::numeric(FLERR, values[4], true, lmp);
// quat[2] = LAMMPS_NS::utils::numeric(FLERR, values[5], true, lmp);
// quat[3] = LAMMPS_NS::utils::numeric(FLERR, values[6], true, lmp);
// LAMMPS_NS::MathExtra::qnormalize(quat);
// reset ellipsoid mass
// previously stored density in rmass
// rmass[m] *= 4.0 * LAMMPS_NS::MY_PI / 3.0 * shape[0] * shape[1] * shape[2];
bonus[nlocal_bonus].ilocal = m;
ellipsoid[m] = nlocal_bonus++;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated bonus memory
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
LAMMPS_NS::bigint AtomVecLM<AtomVecLAMMPS, Dim>::memory_usage_bonus() {
LAMMPS_NS::bigint bytes = 0;
bytes += nmax_bonus * sizeof(Bonus);
return bytes;
}
/* ----------------------------------------------------------------------
initialize non-zero atom quantities
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::create_atom_post(int ilocal) {
rmass[ilocal] = 1.0;
ellipsoid[ilocal] = -1;
}
/* ----------------------------------------------------------------------
modify what AtomVec::data_atom() just unpacked
or initialize other atom quantities
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::data_atom_post(int ilocal) {
ellipsoid_flag = ellipsoid[ilocal];
if (ellipsoid_flag == 0)
ellipsoid_flag = -1;
else if (ellipsoid_flag == 1)
ellipsoid_flag = 0;
else
this->error->one(FLERR,
"Invalid ellipsoid flag in Atoms section of data file");
ellipsoid[ilocal] = ellipsoid_flag;
if (rmass[ilocal] <= 0.0)
this->error->one(FLERR, "Invalid density in Atoms section of data file");
angmom[ilocal][0] = 0.0;
angmom[ilocal][1] = 0.0;
angmom[ilocal][2] = 0.0;
}
/* ----------------------------------------------------------------------
modify values for AtomVec::pack_data() to pack
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::pack_data_pre(int ilocal) {
// double *shape;
ellipsoid_flag = this->atom->ellipsoid[ilocal];
rmass_one = this->atom->rmass[ilocal];
if (ellipsoid_flag < 0)
ellipsoid[ilocal] = 0;
else
ellipsoid[ilocal] = 1;
if (ellipsoid_flag >= 0) {
// shape = bonus[ellipsoid_flag].shape;
// rmass[ilocal] /=
// 4.0 * LAMMPS_NS::MY_PI / 3.0 * shape[0] * shape[1] * shape[2];
}
}
/* ----------------------------------------------------------------------
unmodify values packed by AtomVec::pack_data()
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::pack_data_post(int ilocal) {
ellipsoid[ilocal] = ellipsoid_flag;
rmass[ilocal] = rmass_one;
}
/* ----------------------------------------------------------------------
pack bonus ellipsoid info for writing to data file
if buf is NULL, just return buffer size
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
int AtomVecLM<AtomVecLAMMPS, Dim>::pack_data_bonus(double *, int /*flag*/) {
throw;
// int i, j;
// LAMMPS_NS::tagint *tag = this->atom->tag;
// int nlocal = this->atom->nlocal;
// int m = 0;
// for (i = 0; i < nlocal; i++) {
// if (ellipsoid[i] < 0)
// continue;
// if (buf) {
// buf[m++] = this->ubuf(tag[i]).d;
// j = ellipsoid[i];
// buf[m++] = 2.0 * bonus[j].shape[0];
// buf[m++] = 2.0 * bonus[j].shape[1];
// buf[m++] = 2.0 * bonus[j].shape[2];
// buf[m++] = bonus[j].quat[0];
// buf[m++] = bonus[j].quat[1];
// buf[m++] = bonus[j].quat[2];
// buf[m++] = bonus[j].quat[3];
// } else
// m += this->size_data_bonus;
// }
// return m;
}
/* ----------------------------------------------------------------------
write bonus ellipsoid info to data file
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::write_data_bonus(FILE *, int n, double *,
int /*flag*/) {
int i = 0;
while (i < n) {
// LAMMPS_NS::fmt::print(fp, "{} {} {} {} {} {} {} {}\n",
// this->ubuf(buf[i]).i,
// buf[i + 1], buf[i + 2], buf[i + 3], buf[i + 4],
// buf[i + 5], buf[i + 6], buf[i + 7]);
i += this->size_data_bonus;
}
}
/* ----------------------------------------------------------------------
set shape values in bonus data for particle I
oriented aligned with xyz axes
this may create or delete entry in bonus data
------------------------------------------------------------------------- */
template <typename AtomVecLAMMPS, UInt Dim>
void AtomVecLM<AtomVecLAMMPS, Dim>::set_shape(int i, double shapex,
double shapey, double shapez) {
if (ellipsoid[i] < 0) {
if (shapex == 0.0 && shapey == 0.0 && shapez == 0.0)
return;
if (nlocal_bonus == nmax_bonus)
grow_bonus();
double *shape = bonus[nlocal_bonus].shape;
double *quat = bonus[nlocal_bonus].quat;
shape[0] = shapex;
shape[1] = shapey;
shape[2] = shapez;
quat[0] = 1.0;
quat[1] = 0.0;
quat[2] = 0.0;
quat[3] = 0.0;
bonus[nlocal_bonus].ilocal = i;
ellipsoid[i] = nlocal_bonus++;
} else if (shapex == 0.0 && shapey == 0.0 && shapez == 0.0) {
copy_bonus_all(nlocal_bonus - 1, ellipsoid[i]);
nlocal_bonus--;
ellipsoid[i] = -1;
} else {
double *shape = bonus[ellipsoid[i]].shape;
shape[0] = shapex;
shape[1] = shapey;
shape[2] = shapez;
}
}
__END_LIBMULTISCALE__
#endif //__LIBMULTISCALE_ATOM_VEC_LM_HH__
c4science · Help