Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/filter/compute_centro_symmetry.cc
d8e9d1cc530clammps_2020

Last Change
Raw File

compute_centro_symmetry.cc
View Options

/**
* @file compute_centro_symmetry.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Jul 09 21:59:47 2014
*
* @brief This computes the centro symmetry criterion over a set of points
*
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "compute_centro_symmetry.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "ref_point_data.hh"
#include "trace_atom.hh"
/* -------------------------------------------------------------------------- */
#define SWAP(a, b) \
{ \
tmp = a; \
a = b; \
b = tmp; \
}
#define ISWAP(a, b) \
{ \
itmp = a; \
a = b; \
b = itmp; \
}
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
ComputeCentroSymmetry::ComputeCentroSymmetry(const std::string &name)
: LMObject(name) {
this->createInput("input");
this->createOutput("output");
this->getOutput("output") =
ContainerArray<Real>(this->getID() + ":" + "output");
}
/* -------------------------------------------------------------------------- */
ComputeCentroSymmetry::~ComputeCentroSymmetry() {}
/* -------------------------------------------------------------------------- */
inline void select(UInt k, UInt n, Real *arr) {
UInt i, ir, j, l, mid;
Real a, tmp;
arr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l + 1) {
if (ir == l + 1 && arr[ir] < arr[l]) {
SWAP(arr[l], arr[ir]);
}
return;
} else {
mid = (l + ir) >> 1;
SWAP(arr[mid], arr[l + 1]);
if (arr[l] > arr[ir]) {
SWAP(arr[l], arr[ir]);
}
if (arr[l + 1] > arr[ir]) {
SWAP(arr[l + 1], arr[ir]);
}
if (arr[l] > arr[l + 1]) {
SWAP(arr[l], arr[l + 1]);
}
i = l + 1;
j = ir;
a = arr[l + 1];
for (;;) {
do
i++;
while (arr[i] < a);
do
j--;
while (arr[j] > a);
if (j < i)
break;
SWAP(arr[i], arr[j]);
}
arr[l + 1] = arr[j];
arr[j] = a;
if (j >= k)
ir = j - 1;
if (j <= k)
l = i;
}
}
}
/* ---------------------------------------------------------------------- */
// function which sorts the array arr and move accordingly the iarr entries
//
inline void select2(UInt k, UInt n, Real *arr, UInt *iarr) {
UInt i, ir, j, l, mid, ia, itmp;
Real a, tmp;
arr--;
iarr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l + 1) {
if (ir == l + 1 && arr[ir] < arr[l]) {
SWAP(arr[l], arr[ir]);
ISWAP(iarr[l], iarr[ir]);
}
return;
} else {
mid = (l + ir) >> 1;
SWAP(arr[mid], arr[l + 1]);
ISWAP(iarr[mid], iarr[l + 1]);
if (arr[l] > arr[ir]) {
SWAP(arr[l], arr[ir]);
ISWAP(iarr[l], iarr[ir]);
}
if (arr[l + 1] > arr[ir]) {
SWAP(arr[l + 1], arr[ir]);
ISWAP(iarr[l + 1], iarr[ir]);
}
if (arr[l] > arr[l + 1]) {
SWAP(arr[l], arr[l + 1]);
ISWAP(iarr[l], iarr[l + 1]);
}
i = l + 1;
j = ir;
a = arr[l + 1];
ia = iarr[l + 1];
for (;;) {
do
i++;
while (arr[i] < a);
do
j--;
while (arr[j] > a);
if (j < i)
break;
SWAP(arr[i], arr[j]);
ISWAP(iarr[i], iarr[j]);
}
arr[l + 1] = arr[j];
arr[j] = a;
iarr[l + 1] = iarr[j];
iarr[j] = ia;
if (j >= k)
ir = j - 1;
if (j <= k)
l = i;
}
}
}
/* -------------------------------------------------------------------------- */
template <bool pbc_x, bool pbc_y, bool pbc_z, UInt Dim, typename Cont,
typename ContNeigh>
inline Real ComputeCentroSymmetry::computeCentroSymmetry(
typename Cont::Ref &at, Cont &cont, Real *distsq, UInt *nearest,
UInt maxneigh [[gnu::unused]], Real *pairs, ContNeigh &neighbors) {
Vector<Dim> X = at.position().template cast<Real>();
UInt numneigh = 0;
Cube &bbox = cont.getBoundingBox();
auto d_size = bbox.getSize();
auto d_inv_size_half = bbox.getInvSizeHalf();
UInt tested = 0;
for (auto &at_neighbor : neighbors.to(X)) {
auto X2 = at_neighbor.position();
if (X == X2) {
continue;
}
Real d2 = RefPoint<Dim>::template computeDistance<pbc_x, pbc_y, pbc_z>(
at.position(), X2, bbox);
{
if (CHECK_REF_TRACED(at)) {
DUMP(" testing atom " << at_neighbor, DBG_MESSAGE);
++tested;
}
}
if (d2 > rcut2)
continue; // rejected
DUMP(numneigh << "-th atom selected " << (X[0] - X2[0]) / 3.614999975398666
<< " " << (X[1] - X2[1]) / 3.614999975398666 << " "
<< (X[2] - X2[2]) / 3.614999975398666 << " "
<< " d = " << sqrt(d2) / 3.614999975398666
<< " id = " << at_neighbor.getIndex(),
DBG_DETAIL);
distsq[numneigh] = d2;
nearest[numneigh] = at_neighbor.getIndex();
++numneigh;
}
// if not 12 neighbors, centro = 0.0
if (numneigh < 12)
return 0.0;
LM_ASSERT(numneigh <= maxneigh,
"internal error " << numneigh << " !< " << maxneigh);
DUMP("found " << numneigh << " atoms as first neighbors !", DBG_DETAIL);
// store 12 nearest neighs in 1st 12 locations of distsq and nearest
select2(12, numneigh, distsq, nearest);
{
if (CHECK_REF_TRACED(at)) {
DUMP("CONCERNING atom " << at, DBG_MESSAGE);
DUMP(" I found " << numneigh << "/" << tested << " closest neighbors",
DBG_MESSAGE);
for (UInt j = 0; j < numneigh; j++) {
UInt jj = nearest[j];
typename Cont::Ref atjj = cont.get(jj);
DUMP("d = " << distsq[j] << " " << atjj, DBG_MESSAGE);
}
DUMP("endlist", DBG_MESSAGE);
DUMP("my bbox is " << bbox, DBG_MESSAGE);
}
}
// R = Ri + Rj for each of 66 i,j pairs among 12 neighbors
// pairs = squared length of each R
DUMP("looping over neighbors of atom : " << at, DBG_DETAIL);
typename Cont::Ref atjj;
typename Cont::Ref atkk;
UInt n = 0;
UInt j, jj, k, kk;
for (j = 0; j < 12; j++) {
jj = nearest[j];
atjj = cont.get(jj);
// auto X2 = atjj.position();
for (k = j + 1; k < 12; k++) {
kk = nearest[k];
atkk = cont.get(kk);
pairs[n] = 0.0;
for (UInt i = 0; i < Dim; ++i) {
Real del1 = atjj.position()[i] - X[i];
if ((pbc_x && i == 0) || (pbc_y && i == 1) || (pbc_z && i == 2)) {
int offset = static_cast<int>(del1 * d_inv_size_half[i]);
del1 -= offset * d_size[i];
}
Real del2 = atkk.position()[i] - X[i];
if ((pbc_x && i == 0) || (pbc_y && i == 1) || (pbc_z && i == 2)) {
int offset = static_cast<int>(del2 * d_inv_size_half[i]);
del2 -= offset * d_size[i];
}
Real del = del1 + del2;
pairs[n] += del * del;
}
++n;
}
}
// store 6 smallest pair distances in 1st 6 locations of pairs
select(6, 66, pairs);
// centrosymmetry = sum of 6 smallest squared values
Real value = 0.0;
for (j = 0; j < 6; j++) {
DUMP("value(" << value << ")+= " << pairs[j], DBG_DETAIL);
{
if (CHECK_REF_TRACED(at)) {
DUMP("add " << pairs[j] << " to " << value << " giving "
<< value + pairs[j],
DBG_MESSAGE);
}
}
value += pairs[j];
}
DUMP("centrosymmetry computed for atom " << X << " : " << value, DBG_DETAIL);
return value;
}
/* -------------------------------------------------------------------------- */
template <bool pbc_x, bool pbc_y, bool pbc_z, UInt Dim, typename Cont,
typename ContNeigh>
inline void
ComputeCentroSymmetry::computePerAtomCentroSymmetry(Cont &cont,
ContNeigh &contNeigh) {
UInt index = 0;
UInt maxneigh = 300;
Real *distsq = new Real[maxneigh];
UInt *nearest = new UInt[maxneigh];
Real pairs[66];
for (auto &at : cont) {
if (index % 100000 == 0)
DUMP("computing centrosymmetry for atom " << index << "/" << cont.size(),
DBG_INFO);
Real value = computeCentroSymmetry<pbc_x, pbc_y, pbc_z, Dim>(
at, cont, distsq, nearest, maxneigh, pairs, contNeigh);
this->getOutputAsArray().push_back(value);
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename Cont, typename ContNeigh>
inline void
ComputeCentroSymmetry::computePerAtomCentroSymmetry(Cont &cont,
ContNeigh &contNeigh) {
UInt pbc[3] = {0, 0, 0};
for (UInt i = 0; i < Dim; ++i)
pbc[i] = pbc_flag[i];
if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 0)
computePerAtomCentroSymmetry<false, false, false, Dim>(cont, contNeigh);
if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 0)
computePerAtomCentroSymmetry<true, false, false, Dim>(cont, contNeigh);
if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 0)
computePerAtomCentroSymmetry<true, true, false, Dim>(cont, contNeigh);
if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 1)
computePerAtomCentroSymmetry<true, true, true, Dim>(cont, contNeigh);
if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 1)
computePerAtomCentroSymmetry<true, false, true, Dim>(cont, contNeigh);
if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 0)
computePerAtomCentroSymmetry<false, true, false, Dim>(cont, contNeigh);
if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 1)
computePerAtomCentroSymmetry<false, true, true, Dim>(cont, contNeigh);
if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 1)
computePerAtomCentroSymmetry<false, false, true, Dim>(cont, contNeigh);
}
/* -------------------------------------------------------------------------- */
template <typename Cont> void ComputeCentroSymmetry::build(Cont &cont) {
constexpr UInt Dim = Cont::Dim;
LM_ASSERT(Dim == 3, "cannot work if not in 3D!");
rcut2 = rcut * rcut;
Cube cube = cont.getBoundingBox();
DUMP("compute centro BBOX " << cube, DBG_INFO);
for (UInt i = 0; i < Dim; ++i) {
cube.setXmin(i, cube.getXmin()[i] - this->skin[i]);
cube.setXmax(i, cube.getXmax()[i] + this->skin[i]);
}
this->clear();
using ContNeigh = ContainerNeighborAtoms<typename Cont::Ref, Dim>;
ContNeigh neighbors(this->getID(), cube, rcut);
neighbors.fillGrid(cont);
// neighbors.setPBCFlag(pbc_flag);
#ifdef TRACE_ATOM
UInt indexes[3];
neighbors.coordinates2CartesianIndex(trace_ref, indexes);
DUMP("box indexes containing the traced atom "
<< indexes[0] << " " << indexes[1] << " " << indexes[2],
DBG_MESSAGE);
UInt index;
neighbors.coordinates2Index(trace_ref, index);
DUMP("box index containing the traced atom " << indexes, DBG_MESSAGE);
DUMP("block outside contains " << neighbors.getBlockOutside().size(),
DBG_MESSAGE)
#endif
this->computePerAtomCentroSymmetry<Dim>(cont, neighbors);
}
/* -------------------------------------------------------------------------- */
/* LMDESC CENTRO
This computes the centrosymmetry for a set of points
centro-symmetry criterion is detailed in
Kelchner, C. L., Plimpton, S. J., Hamilton, J. C., Nov 1998. Dislocation
nucleation and defect structure during surface indentation. Phys. Rev. B
58 (17), 11085-11088.
*/
/* LMEXAMPLE
COMPUTE centro CENTRO INPUT md RCUT 3.7
*/
inline void ComputeCentroSymmetry::declareParams() {
ComputeInterface::declareParams();
/* LMKEYWORD RCUT
Specify the cutoff radius to be used for the neighbor detection
*/
this->parseKeyword("RCUT", rcut);
/* LMKEYWORD PBC
Specify the pbc flags (one per direction) for neighbor detection
*/
this->parseVectorKeyword("PBC", spatial_dimension, pbc_flag,
VEC_DEFAULTS(false, false, false));
/* LMKEYWORD SKIN
Specify the additional tolerance to be given per spatial direction to
ensure a large enough search grid
*/
this->parseVectorKeyword("SKIN", spatial_dimension, skin,
VEC_DEFAULTS(0., 0., 0.));
}
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_MAKE_CALL(ComputeCentroSymmetry)
__END_LIBMULTISCALE__
c4science · Help