ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

refactor domain and fix_deform for KOKKOS

Browse Source
This commit is contained in:
alphataubio
2024-08-05 22:01:40 -04:00
parent ab558fb174
commit 88a32e6916
7 changed files with 144 additions and 403 deletions
Download Patch File Download Diff File Expand all files Collapse all files

56
src/KOKKOS/domain_kokkos.cpp
View File

@ -563,7 +563,6 @@ void DomainKokkos::operator()(TagDomain_image_flip, const int &i) const {
void DomainKokkos::lamda2x(int n)
{
atomKK->sync(Device,X_MASK);
x = atomKK->k_x.view<LMPDeviceType>();
copymode = 1;
@ -573,6 +572,21 @@ void DomainKokkos::lamda2x(int n)
atomKK->modified(Device,X_MASK);
}
void DomainKokkos::lamda2x(int n, int groupbit_in)
{
atomKK->sync(Device,X_MASK);
x = atomKK->k_x.view<LMPDeviceType>();
mask = atomKK->k_mask.view<LMPDeviceType>();
mask = atomKK->k_mask.view<LMPDeviceType>();
groupbit = groupbit_in;
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<LMPDeviceType, TagDomain_lamda2x_group>(0,n),*this);
copymode = 0;
atomKK->modified(Device,X_MASK);
}
KOKKOS_INLINE_FUNCTION
void DomainKokkos::operator()(TagDomain_lamda2x, const int &i) const {
const double xi1 = x(i,1);
@ -582,6 +596,17 @@ void DomainKokkos::operator()(TagDomain_lamda2x, const int &i) const {
x(i,2) = h[2]*xi2 + boxlo[2];
}
KOKKOS_INLINE_FUNCTION
void DomainKokkos::operator()(TagDomain_lamda2x_group, const int &i) const {
if (mask[i] & groupbit) {
const double xi1 = x(i,1);
const double xi2 = x(i,2);
x(i,0) = h[0]*x(i,0) + h[5]*xi1 + h[4]*xi2 + boxlo[0];
x(i,1) = h[1]*xi1 + h[3]*xi2 + boxlo[1];
x(i,2) = h[2]*xi2 + boxlo[2];
}
}
/* ----------------------------------------------------------------------
convert box coords to triclinic 0-1 lamda coords for all N atoms
lamda = H^-1 (x - x0)
@ -590,7 +615,6 @@ void DomainKokkos::operator()(TagDomain_lamda2x, const int &i) const {
void DomainKokkos::x2lamda(int n)
{
atomKK->sync(Device,X_MASK);
x = atomKK->k_x.view<LMPDeviceType>();
copymode = 1;
@ -600,6 +624,20 @@ void DomainKokkos::x2lamda(int n)
atomKK->modified(Device,X_MASK);
}
void DomainKokkos::x2lamda(int n, int groupbit_in)
{
atomKK->sync(Device,X_MASK);
x = atomKK->k_x.view<LMPDeviceType>();
mask = atomKK->k_mask.view<LMPDeviceType>();
groupbit = groupbit_in;
copymode = 1;
Kokkos::parallel_for(Kokkos::RangePolicy<LMPDeviceType, TagDomain_x2lamda_group>(0,n),*this);
copymode = 0;
atomKK->modified(Device,X_MASK);
}
KOKKOS_INLINE_FUNCTION
void DomainKokkos::operator()(TagDomain_x2lamda, const int &i) const {
F_FLOAT delta[3];
@ -612,3 +650,17 @@ void DomainKokkos::operator()(TagDomain_x2lamda, const int &i) const {
x(i,2) = h_inv[2]*delta[2];
}
KOKKOS_INLINE_FUNCTION
void DomainKokkos::operator()(TagDomain_x2lamda_group, const int &i) const {
if (mask[i] & groupbit) {
F_FLOAT delta[3];
delta[0] = x(i,0) - boxlo[0];
delta[1] = x(i,1) - boxlo[1];
delta[2] = x(i,2) - boxlo[2];
x(i,0) = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2];
x(i,1) = h_inv[1]*delta[1] + h_inv[3]*delta[2];
x(i,2) = h_inv[2]*delta[2];
}
}

12
src/KOKKOS/domain_kokkos.h
View File

@ -24,7 +24,9 @@ namespace LAMMPS_NS {
struct TagDomain_remap_all{};
struct TagDomain_image_flip{};
struct TagDomain_lamda2x{};
struct TagDomain_lamda2x_group{};
struct TagDomain_x2lamda{};
struct TagDomain_x2lamda_group{};
class DomainKokkos : public Domain {
public:
@ -35,7 +37,9 @@ class DomainKokkos : public Domain {
void remap_all();
void image_flip(int, int, int);
void x2lamda(int) override;
void x2lamda(int,int) override;
void lamda2x(int) override;
void lamda2x(int,int) override;
// forward remaining x2lamda() and lambda2x() variants to parent class
void x2lamda(double *a, double *b) override { Domain::x2lamda(a,b); }
void lamda2x(double *a, double *b) override { Domain::lamda2x(a,b); }
@ -54,18 +58,26 @@ class DomainKokkos : public Domain {
KOKKOS_INLINE_FUNCTION
void operator()(TagDomain_lamda2x, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagDomain_lamda2x_group, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagDomain_x2lamda, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagDomain_x2lamda_group, const int&) const;
static KOKKOS_INLINE_FUNCTION
Few<double,3> unmap(Few<double,3> prd, Few<double,6> h, int triclinic,
Few<double,3> x, imageint image);
private:
int groupbit;
double lo[3],hi[3],period[3];
int n_flip, m_flip, p_flip;
ArrayTypes<LMPDeviceType>::t_x_array x;
ArrayTypes<LMPDeviceType>::t_imageint_1d image;
ArrayTypes<LMPDeviceType>::t_int_1d mask;
};
KOKKOS_INLINE_FUNCTION

367
src/KOKKOS/fix_deform_kokkos.cpp
View File

@ -13,7 +13,8 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Mitch Murphy (alphataubio@gmail.com)
Contributing author: Pieter in 't Veld (SNL)
Refactoring (2024/08): Mitch Murphy (alphataubio@gmail.com)
------------------------------------------------------------------------- */
#include "fix_deform_kokkos.h"
@ -26,7 +27,6 @@
#include "irregular.h"
#include "kspace.h"
#include "math_const.h"
#include "memory_kokkos.h"
#include "modify.h"
#include "update.h"
#include "variable.h"
@ -37,38 +37,16 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{NONE=0,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE,VARIABLE};
enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE};
/* ---------------------------------------------------------------------- */
template <class DeviceType>
FixDeformKokkos<DeviceType>::FixDeformKokkos(LAMMPS *lmp, int narg, char **arg) : FixDeform(lmp, narg, arg)
FixDeformKokkos::FixDeformKokkos(LAMMPS *lmp, int narg, char **arg) : FixDeform(lmp, narg, arg)
{
kokkosable = 1;
atomKK = (AtomKokkos *) atom;
domainKK = (DomainKokkos *) domain;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = EMPTY_MASK;
datamask_modify = EMPTY_MASK;
memoryKK->create_kokkos(k_set,6,"fix_deform:set");
d_set = k_set.template view<DeviceType>();
}
template<class DeviceType>
FixDeformKokkos<DeviceType>::~FixDeformKokkos()
{
if (copymode) return;
memoryKK->destroy_kokkos(k_set,set);
}
template <class DeviceType>
void FixDeformKokkos<DeviceType>::init()
{
FixDeform::init();
memcpy((void *)k_set.h_view.data(), (void *)set, sizeof(Set)*6);
k_set.template modify<LMPHostType>();
k_set.template sync<DeviceType>();
}
/* ----------------------------------------------------------------------
@ -81,335 +59,18 @@ void FixDeformKokkos<DeviceType>::init()
image flags to new values, making eqs in doc of Domain:image_flip incorrect
------------------------------------------------------------------------- */
template <class DeviceType>
void FixDeformKokkos<DeviceType>::pre_exchange()
void FixDeformKokkos::pre_exchange()
{
if (flip == 0) return;
domainKK->yz = d_set[3].tilt_target = d_set[3].tilt_flip;
domainKK->xz = d_set[4].tilt_target = d_set[4].tilt_flip;
domainKK->xy = d_set[5].tilt_target = d_set[5].tilt_flip;
domainKK->set_global_box();
domainKK->set_local_box();
domainKK->image_flip(flipxy, flipxz, flipyz);
// FIXME: just replace with domainKK->remap_all(), is this correct ?
//double **x = atom->x;
//imageint *image = atom->image;
//int nlocal = atom->nlocal;
//for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
//domain->x2lamda(atom->nlocal);
//irregular->migrate_atoms();
//domain->lamda2x(atom->nlocal);
domainKK->remap_all();
flip = 0;
atomKK->sync(Host,ALL_MASK);
FixDeform::pre_exchange();
atomKK->modified(Host,ALL_MASK);
}
/* ---------------------------------------------------------------------- */
template <class DeviceType>
void FixDeformKokkos<DeviceType>::end_of_step()
void FixDeformKokkos::end_of_step()
{
// wrap variable evaluations with clear/add
if (varflag) modify->clearstep_compute();
// set new box size for strain-based dims
apply_strain();
// set new box size for VOLUME dims that are linked to other dims
// NOTE: still need to set h_rate for these dims
apply_volume();
if (varflag) modify->addstep_compute(update->ntimestep + nevery);
update_domain();
// redo KSpace coeffs since box has changed
if (kspace_flag) force->kspace->setup();
atomKK->sync(Host,ALL_MASK);
FixDeform::end_of_step();
atomKK->modified(Host,ALL_MASK);
}
/* ----------------------------------------------------------------------
apply strain controls
------------------------------------------------------------------------- */
template <class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixDeformKokkos<DeviceType>::apply_strain()
{
// for NONE, target is current box size
// for TRATE, set target directly based on current time, also set h_rate
// for WIGGLE, set target directly based on current time, also set h_rate
// for VARIABLE, set target directly via variable eval, also set h_rate
// for others except VOLUME, target is linear value between start and stop
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
for (int i = 0; i < 3; i++) {
if (d_set[i].style == NONE) {
d_set[i].lo_target = domainKK->boxlo[i];
d_set[i].hi_target = domainKK->boxhi[i];
} else if (d_set[i].style == TRATE) {
double delt = (update->ntimestep - update->beginstep) * update->dt;
double shift = 0.5 * ((d_set[i].hi_start - d_set[i].lo_start) * exp(d_set[i].rate * delt));
d_set[i].lo_target = 0.5 * (d_set[i].lo_start + d_set[i].hi_start) - shift;
d_set[i].hi_target = 0.5 * (d_set[i].lo_start + d_set[i].hi_start) + shift;
h_rate[i] = d_set[i].rate * domainKK->h[i];
h_ratelo[i] = -0.5 * h_rate[i];
} else if (d_set[i].style == WIGGLE) {
double delt = (update->ntimestep - update->beginstep) * update->dt;
double shift = 0.5 * d_set[i].amplitude * sin(MY_2PI * delt / d_set[i].tperiod);
d_set[i].lo_target = d_set[i].lo_start - shift;
d_set[i].hi_target = d_set[i].hi_start + shift;
h_rate[i] = MY_2PI / d_set[i].tperiod * d_set[i].amplitude *
cos(MY_2PI * delt / d_set[i].tperiod);
h_ratelo[i] = -0.5 * h_rate[i];
} else if (d_set[i].style == VARIABLE) {
double del = input->variable->compute_equal(d_set[i].hvar);
d_set[i].lo_target = d_set[i].lo_start - 0.5 * del;
d_set[i].hi_target = d_set[i].hi_start + 0.5 * del;
h_rate[i] = input->variable->compute_equal(d_set[i].hratevar);
h_ratelo[i] = -0.5 * h_rate[i];
} else if (d_set[i].style == FINAL || d_set[i].style == DELTA || d_set[i].style == SCALE ||
d_set[i].style == VEL || d_set[i].style == ERATE) {
d_set[i].lo_target = d_set[i].lo_start + delta * (d_set[i].lo_stop - d_set[i].lo_start);
d_set[i].hi_target = d_set[i].hi_start + delta * (d_set[i].hi_stop - d_set[i].hi_start);
}
}
// for triclinic, set new box shape
// for NONE, target is current tilt
// for TRATE, set target directly based on current time. also set h_rate
// for WIGGLE, set target directly based on current time. also set h_rate
// for VARIABLE, set target directly via variable eval. also set h_rate
// for other styles, target is linear value between start and stop values
if (triclinic) {
for (int i = 3; i < 6; i++) {
if (d_set[i].style == NONE) {
if (i == 5) d_set[i].tilt_target = domainKK->xy;
else if (i == 4) d_set[i].tilt_target = domainKK->xz;
else if (i == 3) d_set[i].tilt_target = domainKK->yz;
} else if (d_set[i].style == TRATE) {
double delt = (update->ntimestep - update->beginstep) * update->dt;
d_set[i].tilt_target = d_set[i].tilt_start * exp(d_set[i].rate * delt);
h_rate[i] = d_set[i].rate * domainKK->h[i];
} else if (d_set[i].style == WIGGLE) {
double delt = (update->ntimestep - update->beginstep) * update->dt;
d_set[i].tilt_target = d_set[i].tilt_start +
d_set[i].amplitude * sin(MY_2PI * delt / d_set[i].tperiod);
h_rate[i] = MY_2PI / d_set[i].tperiod * d_set[i].amplitude *
cos(MY_2PI * delt / d_set[i].tperiod);
} else if (d_set[i].style == VARIABLE) {
double delta_tilt = input->variable->compute_equal(d_set[i].hvar);
d_set[i].tilt_target = d_set[i].tilt_start + delta_tilt;
h_rate[i] = input->variable->compute_equal(d_set[i].hratevar);
} else {
d_set[i].tilt_target = d_set[i].tilt_start + delta * (d_set[i].tilt_stop - d_set[i].tilt_start);
}
}
}
}
/* ----------------------------------------------------------------------
apply volume controls
------------------------------------------------------------------------- */
template <class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixDeformKokkos<DeviceType>::apply_volume()
{
for (int i = 0; i < 3; i++) {
if (d_set[i].style != VOLUME) continue;
int dynamic1 = d_set[i].dynamic1;
int dynamic2 = d_set[i].dynamic2;
int fixed = d_set[i].fixed;
double v0 = d_set[i].vol_start;
double shift = 0.0;
if (d_set[i].substyle == ONE_FROM_ONE) {
shift = 0.5 * (v0 / (d_set[dynamic1].hi_target - d_set[dynamic1].lo_target) /
(d_set[fixed].hi_start - d_set[fixed].lo_start));
} else if (d_set[i].substyle == ONE_FROM_TWO) {
shift = 0.5 * (v0 / (d_set[dynamic1].hi_target - d_set[dynamic1].lo_target) /
(d_set[dynamic2].hi_target - d_set[dynamic2].lo_target));
} else if (d_set[i].substyle == TWO_FROM_ONE) {
shift = 0.5 * sqrt(v0 * (d_set[i].hi_start - d_set[i].lo_start) /
(d_set[dynamic1].hi_target - d_set[dynamic1].lo_target) /
(d_set[fixed].hi_start - d_set[fixed].lo_start));
}
h_rate[i] = (2.0 * shift / (domainKK->boxhi[i] - domainKK->boxlo[i]) - 1.0) / update->dt;
h_ratelo[i] = -0.5 * h_rate[i];
d_set[i].lo_target = 0.5 * (d_set[i].lo_start + d_set[i].hi_start) - shift;
d_set[i].hi_target = 0.5 * (d_set[i].lo_start + d_set[i].hi_start) + shift;
}
}
/* ----------------------------------------------------------------------
Update box domain
------------------------------------------------------------------------- */
template <class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixDeformKokkos<DeviceType>::update_domain()
{
// tilt_target can be large positive or large negative value
// add/subtract box lengths until tilt_target is closest to current value
if (triclinic) {
double *h = domainKK->h;
for (int i = 3; i < 6; i++) {
int idenom = 0;
if (i == 5) idenom = 0;
else if (i == 4) idenom = 0;
else if (i == 3) idenom = 1;
double denom = d_set[idenom].hi_target - d_set[idenom].lo_target;
double current = h[i] / h[idenom];
while (d_set[i].tilt_target / denom - current > 0.0)
d_set[i].tilt_target -= denom;
while (d_set[i].tilt_target / denom - current < 0.0)
d_set[i].tilt_target += denom;
if (fabs(d_set[i].tilt_target / denom - 1.0 - current) <
fabs(d_set[i].tilt_target / denom - current))
d_set[i].tilt_target -= denom;
}
}
// if any tilt ratios exceed 0.5, set flip = 1 and compute new tilt values
// do not flip in x or y if non-periodic (can tilt but not flip)
// this is b/c the box length would be changed (dramatically) by flip
// if yz tilt exceeded, adjust C vector by one B vector
// if xz tilt exceeded, adjust C vector by one A vector
// if xy tilt exceeded, adjust B vector by one A vector
// check yz first since it may change xz, then xz check comes after
// flip is performed on next timestep, before reneighboring in pre-exchange()
if (triclinic && flipflag) {
double xprd = d_set[0].hi_target - d_set[0].lo_target;
double yprd = d_set[1].hi_target - d_set[1].lo_target;
double xprdinv = 1.0 / xprd;
double yprdinv = 1.0 / yprd;
if (d_set[3].tilt_target * yprdinv < -0.5 ||
d_set[3].tilt_target * yprdinv > 0.5 ||
d_set[4].tilt_target * xprdinv < -0.5 ||
d_set[4].tilt_target * xprdinv > 0.5 ||
d_set[5].tilt_target * xprdinv < -0.5 ||
d_set[5].tilt_target * xprdinv > 0.5) {
d_set[3].tilt_flip = d_set[3].tilt_target;
d_set[4].tilt_flip = d_set[4].tilt_target;
d_set[5].tilt_flip = d_set[5].tilt_target;
flipxy = flipxz = flipyz = 0;
if (domainKK->yperiodic) {
if (d_set[3].tilt_flip * yprdinv < -0.5) {
d_set[3].tilt_flip += yprd;
d_set[4].tilt_flip += d_set[5].tilt_flip;
flipyz = 1;
} else if (d_set[3].tilt_flip * yprdinv > 0.5) {
d_set[3].tilt_flip -= yprd;
d_set[4].tilt_flip -= d_set[5].tilt_flip;
flipyz = -1;
}
}
if (domainKK->xperiodic) {
if (d_set[4].tilt_flip * xprdinv < -0.5) {
d_set[4].tilt_flip += xprd;
flipxz = 1;
}
if (d_set[4].tilt_flip * xprdinv > 0.5) {
d_set[4].tilt_flip -= xprd;
flipxz = -1;
}
if (d_set[5].tilt_flip * xprdinv < -0.5) {
d_set[5].tilt_flip += xprd;
flipxy = 1;
}
if (d_set[5].tilt_flip * xprdinv > 0.5) {
d_set[5].tilt_flip -= xprd;
flipxy = -1;
}
}
flip = 0;
if (flipxy || flipxz || flipyz) flip = 1;
if (flip) next_reneighbor = update->ntimestep + 1;
}
}
// convert atoms and rigid bodies to lamda coords
if (remapflag == Domain::X_REMAP) {
atomKK->sync(execution_space, X_MASK | MASK_MASK );
d_x = atomKK->k_x.template view<DeviceType>();
d_mask = atomKK->k_mask.template view<DeviceType>();
int nlocal = atomKK->nlocal;
for (int i = 0; i < nlocal; i++)
if (d_mask(i) & groupbit)
domainKK->x2lamda(&d_x(i,0), &d_x(i,0));
for (auto &ifix : rfix)
ifix->deform(0);
}
// reset global and local box to new size/shape
// only if deform fix is controlling the dimension
if (dimflag[0]) {
domainKK->boxlo[0] = d_set[0].lo_target;
domainKK->boxhi[0] = d_set[0].hi_target;
}
if (dimflag[1]) {
domainKK->boxlo[1] = d_set[1].lo_target;
domainKK->boxhi[1] = d_set[1].hi_target;
}
if (dimflag[2]) {
domainKK->boxlo[2] = d_set[2].lo_target;
domainKK->boxhi[2] = d_set[2].hi_target;
}
if (triclinic) {
if (dimflag[3]) domainKK->yz = d_set[3].tilt_target;
if (dimflag[4]) domainKK->xz = d_set[4].tilt_target;
if (dimflag[5]) domainKK->xy = d_set[5].tilt_target;
}
domainKK->set_global_box();
domainKK->set_local_box();
// convert atoms and rigid bodies back to box coords
if (remapflag == Domain::X_REMAP) {
atomKK->sync(execution_space, X_MASK | MASK_MASK );
d_x = atomKK->k_x.template view<DeviceType>();
d_mask = atomKK->k_mask.template view<DeviceType>();
int nlocal = atomKK->nlocal;
for (int i = 0; i < nlocal; i++)
if (d_mask(i) & groupbit)
domainKK->lamda2x(&d_x(i,0), &d_x(i,0));
for (auto &ifix : rfix)
ifix->deform(1);
}
}
namespace LAMMPS_NS {
template class FixDeformKokkos<LMPDeviceType>;
#ifdef LMP_KOKKOS_GPU
template class FixDeformKokkos<LMPHostType>;
#endif
}

32
src/KOKKOS/fix_deform_kokkos.h
View File

@ -13,9 +13,9 @@
#ifdef FIX_CLASS
// clang-format off
FixStyle(deform/kk,FixDeformKokkos<LMPDeviceType>);
FixStyle(deform/kk/device,FixDeformKokkos<LMPDeviceType>);
FixStyle(deform/kk/host,FixDeformKokkos<LMPHostType>);
FixStyle(deform/kk,FixDeformKokkos);
FixStyle(deform/kk/device,FixDeformKokkos);
FixStyle(deform/kk/host,FixDeformKokkos);
// clang-format on
#else
@ -24,40 +24,16 @@ FixStyle(deform/kk/host,FixDeformKokkos<LMPHostType>);
#define LMP_FIX_DEFORM_KOKKOS_H
#include "fix_deform.h"
#include "kokkos_type.h"
namespace LAMMPS_NS {
template <class DeviceType>
class FixDeformKokkos : public FixDeform {
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
FixDeformKokkos(class LAMMPS *, int, char **);
~FixDeformKokkos();
void init() override;
void pre_exchange() override;
void end_of_step() override;
private:
class DomainKokkos *domainKK;
typename AT::t_x_array d_x;
typename AT::t_int_1d d_mask;
Kokkos::DualView<Set*, Kokkos::LayoutRight, LMPDeviceType> k_set;
Kokkos::View<Set*,Kokkos::LayoutRight,DeviceType> d_set;
KOKKOS_INLINE_FUNCTION
void virtual apply_volume();
KOKKOS_INLINE_FUNCTION
void apply_strain();
KOKKOS_INLINE_FUNCTION
void update_domain();
};
}

53
src/domain.cpp
View File

@ -1663,6 +1663,27 @@ void Domain::remap(double *x)
if (triclinic) lamda2x(coord,x);
}
/* ----------------------------------------------------------------------
remap all points into the periodic box no matter how far away
adjust 3 image flags encoded in image accordingly
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
for triclinic, point is converted to lamda coords (0-1) before doing remap
image = 10 bits for each dimension
increment/decrement in wrap-around fashion
------------------------------------------------------------------------- */
void Domain::remap_all()
{
double **x = atom->x;
imageint *image = atom->image;
int nlocal = atom->nlocal;
if (triclinic) x2lamda(nlocal);
for (int i = 0; i < nlocal; i++) remap(x[i],image[i]);
if (triclinic) lamda2x(nlocal);
}
/* ----------------------------------------------------------------------
remap xnew to be within half box length of xold
do it directly, not iteratively, in case is far away
@ -2188,6 +2209,19 @@ void Domain::lamda2x(int n)
}
}
void Domain::lamda2x(int n, int groupbit)
{
double **x = atom->x;
int *mask = atom->mask;
for (int i = 0; i < n; i++)
if (mask[i] & groupbit) {
x[i][0] = h[0]*x[i][0] + h[5]*x[i][1] + h[4]*x[i][2] + boxlo[0];
x[i][1] = h[1]*x[i][1] + h[3]*x[i][2] + boxlo[1];
x[i][2] = h[2]*x[i][2] + boxlo[2];
}
}
/* ----------------------------------------------------------------------
convert box coords to triclinic 0-1 lamda coords for all N atoms
lamda = H^-1 (x - x0)
@ -2209,6 +2243,25 @@ void Domain::x2lamda(int n)
}
}
void Domain::x2lamda(int n, int groupbit)
{
double delta[3];
double **x = atom->x;
int *mask = atom->mask;
for (int i = 0; i < n; i++)
if (mask[i] & groupbit) {
delta[0] = x[i][0] - boxlo[0];
delta[1] = x[i][1] - boxlo[1];
delta[2] = x[i][2] - boxlo[2];
x[i][0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2];
x[i][1] = h_inv[1]*delta[1] + h_inv[3]*delta[2];
x[i][2] = h_inv[2]*delta[2];
}
}
/* ----------------------------------------------------------------------
convert triclinic 0-1 lamda coords to box coords for one atom
x = H lamda + x0;

3
src/domain.h
View File

@ -136,6 +136,7 @@ class Domain : protected Pointers {
void closest_image(const double *const, const double *const, double *const);
void remap(double *, imageint &);
void remap(double *);
void remap_all();
void remap_near(double *, double *);
void unmap_inv(double *x, imageint);
void unmap(double *, imageint);
@ -166,7 +167,9 @@ class Domain : protected Pointers {
void boundary_string(char *);
virtual void lamda2x(int);
virtual void lamda2x(int, int);
virtual void x2lamda(int);
virtual void x2lamda(int, int);
virtual void lamda2x(double *, double *);
virtual void x2lamda(double *, double *);
int inside(double *);

24
src/fix_deform.cpp
View File

@ -648,12 +648,8 @@ void FixDeform::pre_exchange()
domain->set_local_box();
domain->image_flip(flipxy, flipxz, flipyz);
double **x = atom->x;
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
domain->remap_all();
domain->x2lamda(atom->nlocal);
irregular->migrate_atoms();
domain->lamda2x(atom->nlocal);
@ -898,13 +894,7 @@ void FixDeform::update_domain()
// convert atoms and rigid bodies to lamda coords
if (remapflag == Domain::X_REMAP) {
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
domain->x2lamda(x[i], x[i]);
domain->x2lamda(atom->nlocal, groupbit);
for (auto &ifix : rfix)
ifix->deform(0);
@ -937,13 +927,7 @@ void FixDeform::update_domain()
// convert atoms and rigid bodies back to box coords
if (remapflag == Domain::X_REMAP) {
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
domain->lamda2x(x[i], x[i]);
domain->lamda2x(atom->nlocal, groupbit);
for (auto &ifix : rfix)
ifix->deform(1);