renamed files for LAMMPS build system compatibility
This commit is contained in:
pmla
307
src/PTM/compute_ptm_atom.cpp
Normal file
307
src/PTM/compute_ptm_atom.cpp
Normal file
@ -0,0 +1,307 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed
|
||||
under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: PM Larsen (MIT)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <cstdlib>
|
||||
#include <cstring>
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "compute_ptm_atom.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
#include "pair.h"
|
||||
#include "update.h"
|
||||
|
||||
#include "ptm_functions.h"
|
||||
|
||||
#define MAX_NEIGHBORS 30
|
||||
#define NUM_COLUMNS 7
|
||||
#define UNKNOWN 0
|
||||
#define OTHER 8
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
static const char cite_user_ptm_package[] =
|
||||
"USER-PTM package:\n\n"
|
||||
"@Article{larsen2016ptm,\n"
|
||||
" author={Larsen, Peter Mahler and Schmidt, S{\o}ren and Schi{\o}tz, "
|
||||
"Jakob},\n"
|
||||
" title={Robust structural identification via polyhedral template "
|
||||
"matching},\n"
|
||||
" journal={Modelling~Simul.~Mater.~Sci.~Eng.},\n"
|
||||
" year={2016},\n"
|
||||
" number={5},\n"
|
||||
" volume={24},\n"
|
||||
" pages={055007},\n"
|
||||
" DOI = {10.1088/0965-0393/24/5/055007}"
|
||||
"}\n\n";
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputePTMAtom::ComputePTMAtom(LAMMPS *lmp, int narg, char **arg)
|
||||
: Compute(lmp, narg, arg), list(NULL), output(NULL) {
|
||||
if (narg != 5)
|
||||
error->all(FLERR, "Illegal compute ptm/atom command");
|
||||
|
||||
char *structures = arg[3];
|
||||
char *ptr = structures;
|
||||
|
||||
const char *strings[] = {"fcc", "hcp", "bcc", "ico", "sc",
|
||||
"dcub", "dhex", "all", "default"};
|
||||
int32_t flags[] = {
|
||||
PTM_CHECK_FCC,
|
||||
PTM_CHECK_HCP,
|
||||
PTM_CHECK_BCC,
|
||||
PTM_CHECK_ICO,
|
||||
PTM_CHECK_SC,
|
||||
PTM_CHECK_DCUB,
|
||||
PTM_CHECK_DHEX,
|
||||
PTM_CHECK_ALL,
|
||||
PTM_CHECK_FCC | PTM_CHECK_HCP | PTM_CHECK_BCC | PTM_CHECK_ICO};
|
||||
|
||||
input_flags = 0;
|
||||
while (*ptr != '\0') {
|
||||
|
||||
bool found = false;
|
||||
for (int i = 0; i < 9; i++) {
|
||||
int len = strlen(strings[i]);
|
||||
if (strncmp(ptr, strings[i], len) == 0) {
|
||||
input_flags |= flags[i];
|
||||
ptr += len;
|
||||
found = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (!found)
|
||||
error->all(FLERR,
|
||||
"Illegal compute ptm/atom command (invalid structure type)");
|
||||
|
||||
if (*ptr == '\0')
|
||||
break;
|
||||
|
||||
if (*ptr != '-')
|
||||
error->all(FLERR,
|
||||
"Illegal compute ptm/atom command (invalid structure type)");
|
||||
|
||||
ptr++;
|
||||
}
|
||||
|
||||
double threshold = force->numeric(FLERR, arg[4]);
|
||||
if (threshold < 0.0)
|
||||
error->all(FLERR,
|
||||
"Illegal compute ptm/atom command (threshold is negative)");
|
||||
rmsd_threshold = threshold;
|
||||
if (rmsd_threshold == 0)
|
||||
rmsd_threshold = INFINITY;
|
||||
|
||||
peratom_flag = 1;
|
||||
size_peratom_cols = NUM_COLUMNS;
|
||||
create_attribute = 1;
|
||||
nmax = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputePTMAtom::~ComputePTMAtom() { memory->destroy(output); }
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputePTMAtom::init() {
|
||||
if (force->pair == NULL)
|
||||
error->all(FLERR, "Compute ptm/atom requires a pair style be defined");
|
||||
|
||||
int count = 0;
|
||||
for (int i = 0; i < modify->ncompute; i++)
|
||||
if (strcmp(modify->compute[i]->style, "ptm/atom") == 0)
|
||||
count++;
|
||||
if (count > 1 && comm->me == 0)
|
||||
error->warning(FLERR, "More than one compute ptm/atom defined");
|
||||
|
||||
// need an occasional full neighbor list
|
||||
|
||||
int irequest = neighbor->request(this, instance_me);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->compute = 1;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
neighbor->requests[irequest]->occasional = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputePTMAtom::init_list(int id, NeighList *ptr) { list = ptr; }
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
typedef struct {
|
||||
int index;
|
||||
double d;
|
||||
} ptmnbr_t;
|
||||
|
||||
static bool sorthelper_compare(ptmnbr_t const &a, ptmnbr_t const &b) {
|
||||
return a.d < b.d;
|
||||
}
|
||||
|
||||
static int get_neighbors(double *pos, int jnum, int *jlist, double **x,
|
||||
double (*nbr)[3]) {
|
||||
|
||||
ptmnbr_t *nbr_order = new ptmnbr_t[jnum];
|
||||
|
||||
for (int jj = 0; jj < jnum; jj++) {
|
||||
int j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
double dx = pos[0] - x[j][0];
|
||||
double dy = pos[1] - x[j][1];
|
||||
double dz = pos[2] - x[j][2];
|
||||
double rsq = dx * dx + dy * dy + dz * dz;
|
||||
|
||||
nbr_order[jj].index = j;
|
||||
nbr_order[jj].d = rsq;
|
||||
}
|
||||
|
||||
std::sort(nbr_order, nbr_order + jnum, &sorthelper_compare);
|
||||
int num_nbrs = std::min(MAX_NEIGHBORS, jnum);
|
||||
|
||||
nbr[0][0] = nbr[0][1] = nbr[0][2] = 0;
|
||||
for (int jj = 0; jj < num_nbrs; jj++) {
|
||||
|
||||
int j = nbr_order[jj].index;
|
||||
nbr[jj + 1][0] = x[j][0] - pos[0];
|
||||
nbr[jj + 1][1] = x[j][1] - pos[1];
|
||||
nbr[jj + 1][2] = x[j][2] - pos[2];
|
||||
}
|
||||
|
||||
delete[] nbr_order;
|
||||
return num_nbrs;
|
||||
}
|
||||
|
||||
void ComputePTMAtom::compute_peratom() {
|
||||
// PTM global initialization. If already initialized this function does
|
||||
// nothing.
|
||||
ptm_initialize_global();
|
||||
|
||||
// initialize PTM local storage
|
||||
ptm_local_handle_t local_handle = ptm_initialize_local();
|
||||
|
||||
invoked_peratom = update->ntimestep;
|
||||
|
||||
// grow arrays if necessary
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(output);
|
||||
nmax = atom->nmax;
|
||||
|
||||
memory->create(output, nmax, NUM_COLUMNS, "ptm:ptm_output");
|
||||
array_atom = output;
|
||||
}
|
||||
|
||||
// invoke full neighbor list (will copy or build if necessary)
|
||||
neighbor->build_one(list);
|
||||
|
||||
int inum = list->inum;
|
||||
int *ilist = list->ilist;
|
||||
int *numneigh = list->numneigh;
|
||||
int **firstneigh = list->firstneigh;
|
||||
|
||||
double **x = atom->x;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int ii = 0; ii < inum; ii++) {
|
||||
|
||||
int i = ilist[ii];
|
||||
output[i][0] = UNKNOWN;
|
||||
if (!(mask[i] & groupbit))
|
||||
continue;
|
||||
|
||||
double *pos = x[i];
|
||||
|
||||
int *jlist = firstneigh[i];
|
||||
int jnum = numneigh[i];
|
||||
if (jnum <= 0)
|
||||
continue;
|
||||
|
||||
// get neighbours ordered by increasing distance
|
||||
double nbr[MAX_NEIGHBORS + 1][3];
|
||||
int num_nbrs = get_neighbors(pos, jnum, jlist, x, nbr);
|
||||
|
||||
// check that we have enough neighbours for the desired structure types
|
||||
int32_t flags = 0;
|
||||
if (num_nbrs >= PTM_NUM_NBRS_SC && (input_flags & PTM_CHECK_SC))
|
||||
flags |= PTM_CHECK_SC;
|
||||
if (num_nbrs >= PTM_NUM_NBRS_FCC && (input_flags & PTM_CHECK_FCC))
|
||||
flags |= PTM_CHECK_FCC;
|
||||
if (num_nbrs >= PTM_NUM_NBRS_HCP && (input_flags & PTM_CHECK_HCP))
|
||||
flags |= PTM_CHECK_HCP;
|
||||
if (num_nbrs >= PTM_NUM_NBRS_ICO && (input_flags & PTM_CHECK_ICO))
|
||||
flags |= PTM_CHECK_ICO;
|
||||
if (num_nbrs >= PTM_NUM_NBRS_BCC && (input_flags & PTM_CHECK_BCC))
|
||||
flags |= PTM_CHECK_BCC;
|
||||
if (num_nbrs >= PTM_NUM_NBRS_DCUB && (input_flags & PTM_CHECK_DCUB))
|
||||
flags |= PTM_CHECK_DCUB;
|
||||
if (num_nbrs >= PTM_NUM_NBRS_DHEX && (input_flags & PTM_CHECK_DHEX))
|
||||
flags |= PTM_CHECK_DHEX;
|
||||
|
||||
// now run PTM
|
||||
int8_t mapping[MAX_NEIGHBORS + 1];
|
||||
int32_t type, alloy_type;
|
||||
double scale, rmsd, interatomic_distance, lattice_constant;
|
||||
double q[4], F[9], F_res[3], U[9], P[9];
|
||||
ptm_index(local_handle, flags, num_nbrs + 1, nbr, NULL, true, &type,
|
||||
&alloy_type, &scale, &rmsd, q, F, F_res, U, P, mapping,
|
||||
&interatomic_distance, &lattice_constant);
|
||||
|
||||
if (rmsd > rmsd_threshold) {
|
||||
type = PTM_MATCH_NONE;
|
||||
}
|
||||
|
||||
// printf("%d type=%d rmsd=%f\n", i, type, rmsd);
|
||||
|
||||
if (type == PTM_MATCH_NONE)
|
||||
type = OTHER;
|
||||
|
||||
output[i][0] = type;
|
||||
output[i][1] = rmsd;
|
||||
output[i][2] = interatomic_distance;
|
||||
output[i][3] = q[0];
|
||||
output[i][4] = q[1];
|
||||
output[i][5] = q[2];
|
||||
output[i][6] = q[3];
|
||||
}
|
||||
|
||||
// printf("finished ptm analysis\n");
|
||||
ptm_uninitialize_local(local_handle);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
memory usage of local atom-based array
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double ComputePTMAtom::memory_usage() {
|
||||
double bytes = nmax * NUM_COLUMNS * sizeof(double);
|
||||
bytes += nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
48
src/PTM/compute_ptm_atom.h
Normal file
48
src/PTM/compute_ptm_atom.h
Normal file
@ -0,0 +1,48 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
|
||||
ComputeStyle(ptm/atom,ComputePTMAtom)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_PTM_ATOM_H
|
||||
#define LMP_COMPUTE_PTM_ATOM_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputePTMAtom : public Compute {
|
||||
public:
|
||||
ComputePTMAtom(class LAMMPS *, int, char **);
|
||||
~ComputePTMAtom();
|
||||
void init();
|
||||
void init_list(int, class NeighList *);
|
||||
void compute_peratom();
|
||||
double memory_usage();
|
||||
|
||||
private:
|
||||
int nmax;
|
||||
int32_t input_flags;
|
||||
double rmsd_threshold;
|
||||
class NeighList *list;
|
||||
double **output;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
101
src/PTM/ptm_alloy_types.cpp
Normal file
101
src/PTM/ptm_alloy_types.cpp
Normal file
@ -0,0 +1,101 @@
|
||||
#include <algorithm>
|
||||
#include "ptm_constants.h"
|
||||
#include "ptm_initialize_data.h"
|
||||
|
||||
|
||||
#define NUM_ALLOY_TYPES 3
|
||||
static uint32_t typedata[NUM_ALLOY_TYPES][3] = {
|
||||
{PTM_MATCH_FCC, PTM_ALLOY_L10, 0x000001fe},
|
||||
{PTM_MATCH_FCC, PTM_ALLOY_L12_CU, 0x0000001e},
|
||||
{PTM_MATCH_FCC, PTM_ALLOY_L12_AU, 0x00001ffe},
|
||||
};
|
||||
|
||||
static bool test_pure(int num_nbrs, int32_t* numbers)
|
||||
{
|
||||
for (int i=1;i<num_nbrs + 1;i++)
|
||||
if (numbers[i] != numbers[0])
|
||||
return false;
|
||||
return true;
|
||||
}
|
||||
|
||||
static bool test_binary(int num_nbrs, int32_t* numbers)
|
||||
{
|
||||
int a = numbers[0], b = -1;
|
||||
for (int i=1;i<num_nbrs + 1;i++)
|
||||
{
|
||||
if (numbers[i] != a)
|
||||
{
|
||||
if (b == -1)
|
||||
b = numbers[i];
|
||||
else if (numbers[i] != b)
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
static bool test_shell_structure(const refdata_t* ref, int8_t* mapping, int32_t* numbers, int num_inner)
|
||||
{
|
||||
int8_t binary[PTM_MAX_POINTS];
|
||||
for (int i=0;i<ref->num_nbrs+1;i++)
|
||||
binary[i] = numbers[mapping[i]] == numbers[0] ? 0 : 1;
|
||||
|
||||
for (int i=1;i<num_inner + 1;i++)
|
||||
if (binary[i] == binary[0])
|
||||
return false;
|
||||
|
||||
for (int i=num_inner+1;i<ref->num_nbrs+1;i++)
|
||||
if (binary[i] != binary[0])
|
||||
return false;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
static int32_t canonical_alloy_representation(const refdata_t* ref, int8_t* mapping, int32_t* numbers)
|
||||
{
|
||||
int8_t binary[PTM_MAX_POINTS];
|
||||
for (int i=0;i<ref->num_nbrs+1;i++)
|
||||
binary[i] = numbers[mapping[i]] == numbers[0] ? 0 : 1;
|
||||
|
||||
int8_t temp[PTM_MAX_POINTS];
|
||||
uint32_t best = 0xFFFFFFFF;
|
||||
for (int j=0;j<ref->num_mappings;j++)
|
||||
{
|
||||
for (int i=0;i<ref->num_nbrs+1;i++)
|
||||
temp[ref->mapping[j][i]] = binary[i];
|
||||
|
||||
uint32_t code = 0;
|
||||
for (int i=0;i<ref->num_nbrs+1;i++)
|
||||
code |= (temp[i] << i);
|
||||
|
||||
best = std::min(best, code);
|
||||
}
|
||||
|
||||
return best;
|
||||
}
|
||||
|
||||
int32_t find_alloy_type(const refdata_t* ref, int8_t* mapping, int32_t* numbers)
|
||||
{
|
||||
if (test_pure(ref->num_nbrs, numbers))
|
||||
return PTM_ALLOY_PURE;
|
||||
|
||||
if (!test_binary(ref->num_nbrs, numbers))
|
||||
return PTM_ALLOY_NONE;
|
||||
|
||||
uint32_t code = canonical_alloy_representation(ref, mapping, numbers);
|
||||
for (int i=0;i<NUM_ALLOY_TYPES;i++)
|
||||
if ((uint32_t)ref->type == typedata[i][0] && code == typedata[i][2])
|
||||
return typedata[i][1];
|
||||
|
||||
if (ref->type == PTM_MATCH_BCC)
|
||||
if (test_shell_structure(ref, mapping, numbers, 8))
|
||||
return PTM_ALLOY_B2;
|
||||
|
||||
if (ref->type == PTM_MATCH_DCUB || ref->type == PTM_MATCH_DHEX)
|
||||
if (test_shell_structure(ref, mapping, numbers, 4))
|
||||
return PTM_ALLOY_SIC;
|
||||
|
||||
return PTM_ALLOY_NONE;
|
||||
}
|
||||
|
||||
9
src/PTM/ptm_alloy_types.h
Normal file
9
src/PTM/ptm_alloy_types.h
Normal file
@ -0,0 +1,9 @@
|
||||
#ifndef PTM_ALLOY_TYPES_H
|
||||
#define PTM_ALLOY_TYPES_H
|
||||
|
||||
#include "ptm_initialize_data.h"
|
||||
|
||||
int32_t find_alloy_type(const refdata_t* ref, int8_t* mapping, int32_t* numbers);
|
||||
|
||||
#endif
|
||||
|
||||
167
src/PTM/ptm_canonical_coloured.cpp
Normal file
167
src/PTM/ptm_canonical_coloured.cpp
Normal file
@ -0,0 +1,167 @@
|
||||
#include <string.h>
|
||||
#include <climits>
|
||||
#include <algorithm>
|
||||
#include "ptm_graph_tools.h"
|
||||
#include "ptm_constants.h"
|
||||
|
||||
|
||||
static bool weinberg_coloured(int num_nodes, int num_edges, int8_t common[PTM_MAX_NBRS][PTM_MAX_NBRS], int8_t* colours, int8_t* best_code, int8_t* canonical_labelling, int a, int b)
|
||||
{
|
||||
bool m[PTM_MAX_NBRS][PTM_MAX_NBRS];
|
||||
memset(m, 0, sizeof(bool) * PTM_MAX_NBRS * PTM_MAX_NBRS);
|
||||
|
||||
int8_t index[PTM_MAX_NBRS];
|
||||
memset(index, -1, sizeof(int8_t) * PTM_MAX_NBRS);
|
||||
|
||||
|
||||
int n = 0;
|
||||
index[a] = colours[a] * num_nodes + n++;
|
||||
if (index[a] > best_code[0])
|
||||
return false;
|
||||
|
||||
bool winning = false;
|
||||
if (index[a] < best_code[0])
|
||||
{
|
||||
best_code[0] = index[a];
|
||||
winning = true;
|
||||
}
|
||||
|
||||
int c = -1;
|
||||
for (int it=1;it<2*num_edges;it++)
|
||||
{
|
||||
bool newvertex = index[b] == -1;
|
||||
|
||||
if (newvertex)
|
||||
index[b] = colours[b] * num_nodes + n++;
|
||||
|
||||
if (!winning && index[b] > best_code[it])
|
||||
return false;
|
||||
|
||||
if (winning || index[b] < best_code[it])
|
||||
{
|
||||
winning = true;
|
||||
best_code[it] = index[b];
|
||||
}
|
||||
|
||||
if (newvertex)
|
||||
{
|
||||
//When a new vertex is reached, take the right-most edge
|
||||
//relative to the edge on which the vertex is reached.
|
||||
|
||||
c = common[a][b];
|
||||
}
|
||||
else if (m[b][a] == false)
|
||||
{
|
||||
//When an old vertex is reached on a new path, go back
|
||||
//in the opposite direction.
|
||||
|
||||
c = a;
|
||||
}
|
||||
else
|
||||
{
|
||||
//When an old vertex is reached on an old path, leave the
|
||||
//vertex on the right-most edge that has not previously
|
||||
//been traversed in that direction.
|
||||
|
||||
c = common[a][b];
|
||||
while (m[b][c] == true)
|
||||
c = common[c][b];
|
||||
}
|
||||
|
||||
m[a][b] = true;
|
||||
a = b;
|
||||
b = c;
|
||||
}
|
||||
|
||||
if (winning)
|
||||
{
|
||||
memcpy(canonical_labelling, index, sizeof(int8_t) * num_nodes);
|
||||
return true;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
int canonical_form_coloured(int num_facets, int8_t facets[][3], int num_nodes, int8_t* degree, int8_t* colours, int8_t* canonical_labelling, int8_t* best_code, uint64_t* p_hash)
|
||||
{
|
||||
int8_t common[PTM_MAX_NBRS][PTM_MAX_NBRS] = {{0}};
|
||||
int num_edges = 3 * num_facets / 2;
|
||||
if (!build_facet_map(num_facets, facets, common))
|
||||
return -1;
|
||||
|
||||
memset(best_code, SCHAR_MAX, sizeof(int8_t) * 2 * PTM_MAX_EDGES);
|
||||
|
||||
bool equal = true;
|
||||
for (int i = 1;i<num_nodes;i++)
|
||||
if (degree[i] != degree[0] || colours[i] != colours[0])
|
||||
equal = false;
|
||||
|
||||
if (equal)
|
||||
{
|
||||
weinberg_coloured(num_nodes, num_edges, common, colours, best_code, canonical_labelling, facets[0][0], facets[0][1]);
|
||||
}
|
||||
else
|
||||
{
|
||||
uint32_t best_degree = 0;
|
||||
for (int i = 0;i<num_facets;i++)
|
||||
{
|
||||
int a = facets[i][0];
|
||||
int b = facets[i][1];
|
||||
int c = facets[i][2];
|
||||
|
||||
//int da = colours[a] * num_nodes + degree[a];
|
||||
//int db = colours[b] * num_nodes + degree[b];
|
||||
//int dc = colours[c] * num_nodes + degree[c];
|
||||
|
||||
int da = degree[a];
|
||||
int db = degree[b];
|
||||
int dc = degree[c];
|
||||
|
||||
best_degree = std::max(best_degree, ((uint32_t)da << 16) | ((uint32_t)db << 8) | ((uint32_t)dc << 0));
|
||||
best_degree = std::max(best_degree, ((uint32_t)da << 0) | ((uint32_t)db << 16) | ((uint32_t)dc << 8));
|
||||
best_degree = std::max(best_degree, ((uint32_t)da << 8) | ((uint32_t)db << 0) | ((uint32_t)dc << 16));
|
||||
}
|
||||
|
||||
for (int i = 0;i<num_facets;i++)
|
||||
{
|
||||
int a = facets[i][0];
|
||||
int b = facets[i][1];
|
||||
int c = facets[i][2];
|
||||
|
||||
//int da = colours[a] * num_nodes + degree[a];
|
||||
//int db = colours[b] * num_nodes + degree[b];
|
||||
//int dc = colours[c] * num_nodes + degree[c];
|
||||
|
||||
int da = degree[a];
|
||||
int db = degree[b];
|
||||
int dc = degree[c];
|
||||
|
||||
if (best_degree == (((uint32_t)da << 16) | ((uint32_t)db << 8) | ((uint32_t)dc << 0)))
|
||||
weinberg_coloured(num_nodes, num_edges, common, colours, best_code, canonical_labelling, a, b);
|
||||
|
||||
if (best_degree == (((uint32_t)da << 0) | ((uint32_t)db << 16) | ((uint32_t)dc << 8)))
|
||||
weinberg_coloured(num_nodes, num_edges, common, colours, best_code, canonical_labelling, b, c);
|
||||
|
||||
if (best_degree == (((uint32_t)da << 8) | ((uint32_t)db << 0) | ((uint32_t)dc << 16)))
|
||||
weinberg_coloured(num_nodes, num_edges, common, colours, best_code, canonical_labelling, c, a);
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = num_nodes-1;i>=0;i--)
|
||||
canonical_labelling[i+1] = (canonical_labelling[i] % num_nodes) + 1;
|
||||
canonical_labelling[0] = 0;
|
||||
|
||||
uint64_t hash = 0;
|
||||
for (int i = 0;i<2 * num_edges;i++)
|
||||
{
|
||||
uint64_t e = best_code[i];
|
||||
e += i % 8;
|
||||
e &= 0xF;
|
||||
e <<= (4 * i) % 64;
|
||||
hash ^= e;
|
||||
}
|
||||
|
||||
*p_hash = hash;
|
||||
return PTM_NO_ERROR;
|
||||
}
|
||||
|
||||
9
src/PTM/ptm_canonical_coloured.h
Normal file
9
src/PTM/ptm_canonical_coloured.h
Normal file
@ -0,0 +1,9 @@
|
||||
#ifndef PTM_CANONICAL_COLOURED_H
|
||||
#define PTM_CANONICAL_COLOURED_H
|
||||
|
||||
#include <stdint.h>
|
||||
|
||||
int canonical_form_coloured(int num_facets, int8_t facets[][3], int num_nodes, int8_t* degree, int8_t* colours, int8_t* canonical_labelling, int8_t* best_code, uint64_t* p_hash);
|
||||
|
||||
#endif
|
||||
|
||||
174
src/PTM/ptm_constants.h
Normal file
174
src/PTM/ptm_constants.h
Normal file
@ -0,0 +1,174 @@
|
||||
#ifndef PTM_CONSTANTS_H
|
||||
#define PTM_CONSTANTS_H
|
||||
|
||||
//------------------------------------
|
||||
// definitions
|
||||
//------------------------------------
|
||||
#define PTM_NO_ERROR 0
|
||||
|
||||
|
||||
#define PTM_CHECK_FCC (1 << 0)
|
||||
#define PTM_CHECK_HCP (1 << 1)
|
||||
#define PTM_CHECK_BCC (1 << 2)
|
||||
#define PTM_CHECK_ICO (1 << 3)
|
||||
#define PTM_CHECK_SC (1 << 4)
|
||||
#define PTM_CHECK_DCUB (1 << 5)
|
||||
#define PTM_CHECK_DHEX (1 << 6)
|
||||
#define PTM_CHECK_NONDIAMOND (PTM_CHECK_SC | PTM_CHECK_FCC | PTM_CHECK_HCP | PTM_CHECK_ICO | PTM_CHECK_BCC)
|
||||
#define PTM_CHECK_ALL (PTM_CHECK_SC | PTM_CHECK_FCC | PTM_CHECK_HCP | PTM_CHECK_ICO | PTM_CHECK_BCC | PTM_CHECK_DCUB | PTM_CHECK_DHEX)
|
||||
|
||||
#define PTM_MATCH_NONE 0
|
||||
#define PTM_MATCH_FCC 1
|
||||
#define PTM_MATCH_HCP 2
|
||||
#define PTM_MATCH_BCC 3
|
||||
#define PTM_MATCH_ICO 4
|
||||
#define PTM_MATCH_SC 5
|
||||
#define PTM_MATCH_DCUB 6
|
||||
#define PTM_MATCH_DHEX 7
|
||||
|
||||
#define PTM_ALLOY_NONE 0
|
||||
#define PTM_ALLOY_PURE 1
|
||||
#define PTM_ALLOY_L10 2
|
||||
#define PTM_ALLOY_L12_CU 3
|
||||
#define PTM_ALLOY_L12_AU 4
|
||||
#define PTM_ALLOY_B2 5
|
||||
#define PTM_ALLOY_SIC 6
|
||||
|
||||
|
||||
#define PTM_MAX_INPUT_POINTS 35
|
||||
#define PTM_MAX_NBRS 16
|
||||
#define PTM_MAX_POINTS (PTM_MAX_NBRS + 1)
|
||||
#define PTM_MAX_FACETS 28 //2 * PTM_MAX_NBRS - 4
|
||||
#define PTM_MAX_EDGES 42 //3 * PTM_MAX_NBRS - 6
|
||||
|
||||
|
||||
//------------------------------------
|
||||
// number of neighbours
|
||||
//------------------------------------
|
||||
#define PTM_NUM_NBRS_FCC 12
|
||||
#define PTM_NUM_NBRS_HCP 12
|
||||
#define PTM_NUM_NBRS_BCC 14
|
||||
#define PTM_NUM_NBRS_ICO 12
|
||||
#define PTM_NUM_NBRS_SC 6
|
||||
#define PTM_NUM_NBRS_DCUB 16
|
||||
#define PTM_NUM_NBRS_DHEX 16
|
||||
|
||||
#define PTM_NUM_POINTS_FCC (PTM_NUM_NBRS_FCC + 1)
|
||||
#define PTM_NUM_POINTS_HCP (PTM_NUM_NBRS_HCP + 1)
|
||||
#define PTM_NUM_POINTS_BCC (PTM_NUM_NBRS_BCC + 1)
|
||||
#define PTM_NUM_POINTS_ICO (PTM_NUM_NBRS_ICO + 1)
|
||||
#define PTM_NUM_POINTS_SC (PTM_NUM_NBRS_SC + 1)
|
||||
#define PTM_NUM_POINTS_DCUB (PTM_NUM_NBRS_DCUB + 1)
|
||||
#define PTM_NUM_POINTS_DHEX (PTM_NUM_NBRS_DHEX + 1)
|
||||
|
||||
const int ptm_num_nbrs[8] = {0, PTM_NUM_NBRS_FCC, PTM_NUM_NBRS_HCP, PTM_NUM_NBRS_BCC, PTM_NUM_NBRS_ICO, PTM_NUM_NBRS_SC, PTM_NUM_NBRS_DCUB, PTM_NUM_NBRS_DHEX};
|
||||
|
||||
//------------------------------------
|
||||
// template structures
|
||||
//------------------------------------
|
||||
|
||||
//these point sets have barycentre {0, 0, 0} and are scaled such that the mean neighbour distance is 1
|
||||
|
||||
const double ptm_template_fcc[PTM_NUM_POINTS_FCC][3] = { { 0. , 0. , 0. },
|
||||
{ 0. , 0.707106781187, 0.707106781187 },
|
||||
{ 0. , -0.707106781187, -0.707106781187 },
|
||||
{ 0. , 0.707106781187, -0.707106781187 },
|
||||
{ 0. , -0.707106781187, 0.707106781187 },
|
||||
{ 0.707106781187, 0. , 0.707106781187 },
|
||||
{ -0.707106781187, 0. , -0.707106781187 },
|
||||
{ 0.707106781187, 0. , -0.707106781187 },
|
||||
{ -0.707106781187, 0. , 0.707106781187 },
|
||||
{ 0.707106781187, 0.707106781187, 0. },
|
||||
{ -0.707106781187, -0.707106781187, 0. },
|
||||
{ 0.707106781187, -0.707106781187, 0. },
|
||||
{ -0.707106781187, 0.707106781187, 0. } };
|
||||
|
||||
const double ptm_template_hcp[PTM_NUM_POINTS_HCP][3] = { { 0. , 0. , 0. },
|
||||
{ 0.707106781186, 0. , 0.707106781186 },
|
||||
{ -0.235702260395, -0.942809041583, -0.235702260395 },
|
||||
{ 0.707106781186, 0.707106781186, 0. },
|
||||
{ -0.235702260395, -0.235702260395, -0.942809041583 },
|
||||
{ 0. , 0.707106781186, 0.707106781186 },
|
||||
{ -0.942809041583, -0.235702260395, -0.235702260395 },
|
||||
{ -0.707106781186, 0.707106781186, 0. },
|
||||
{ 0. , 0.707106781186, -0.707106781186 },
|
||||
{ 0.707106781186, 0. , -0.707106781186 },
|
||||
{ 0.707106781186, -0.707106781186, 0. },
|
||||
{ -0.707106781186, 0. , 0.707106781186 },
|
||||
{ 0. , -0.707106781186, 0.707106781186 } };
|
||||
|
||||
const double ptm_template_bcc[PTM_NUM_POINTS_BCC][3] = { { 0. , 0. , 0. },
|
||||
{ -0.541451884327, -0.541451884327, -0.541451884327 },
|
||||
{ 0.541451884327, 0.541451884327, 0.541451884327 },
|
||||
{ 0.541451884327, -0.541451884327, -0.541451884327 },
|
||||
{ -0.541451884327, 0.541451884327, 0.541451884327 },
|
||||
{ -0.541451884327, 0.541451884327, -0.541451884327 },
|
||||
{ 0.541451884327, -0.541451884327, 0.541451884327 },
|
||||
{ -0.541451884327, -0.541451884327, 0.541451884327 },
|
||||
{ 0.541451884327, 0.541451884327, -0.541451884327 },
|
||||
{ 0. , 0. , -1.082903768655 },
|
||||
{ 0. , 0. , 1.082903768655 },
|
||||
{ 0. , -1.082903768655, 0. },
|
||||
{ 0. , 1.082903768655, 0. },
|
||||
{ -1.082903768655, 0. , 0. },
|
||||
{ 1.082903768655, 0. , 0. } };
|
||||
|
||||
const double ptm_template_ico[PTM_NUM_POINTS_ICO][3] = { { 0. , 0. , 0. },
|
||||
{ 0. , 0.525731112119, 0.850650808352 },
|
||||
{ 0. , -0.525731112119, -0.850650808352 },
|
||||
{ 0. , 0.525731112119, -0.850650808352 },
|
||||
{ 0. , -0.525731112119, 0.850650808352 },
|
||||
{ -0.525731112119, -0.850650808352, 0. },
|
||||
{ 0.525731112119, 0.850650808352, 0. },
|
||||
{ 0.525731112119, -0.850650808352, 0. },
|
||||
{ -0.525731112119, 0.850650808352, 0. },
|
||||
{ -0.850650808352, 0. , -0.525731112119 },
|
||||
{ 0.850650808352, 0. , 0.525731112119 },
|
||||
{ 0.850650808352, 0. , -0.525731112119 },
|
||||
{ -0.850650808352, 0. , 0.525731112119 } };
|
||||
|
||||
const double ptm_template_sc[PTM_NUM_POINTS_SC][3] = { { 0. , 0. , 0. },
|
||||
{ 0. , 0. , -1. },
|
||||
{ 0. , 0. , 1. },
|
||||
{ 0. , -1. , 0. },
|
||||
{ 0. , 1. , 0. },
|
||||
{ -1. , 0. , 0. },
|
||||
{ 1. , 0. , 0. } };
|
||||
|
||||
const double ptm_template_dcub[PTM_NUM_POINTS_DCUB][3] = { { 0. , 0. , 0. },
|
||||
{ -0.391491627053, 0.391491627053, 0.391491627053 },
|
||||
{ -0.391491627053, -0.391491627053, -0.391491627053 },
|
||||
{ 0.391491627053, -0.391491627053, 0.391491627053 },
|
||||
{ 0.391491627053, 0.391491627053, -0.391491627053 },
|
||||
{ -0.782983254107, 0. , 0.782983254107 },
|
||||
{ -0.782983254107, 0.782983254107, 0. },
|
||||
{ 0. , 0.782983254107, 0.782983254107 },
|
||||
{ -0.782983254107, -0.782983254107, 0. },
|
||||
{ -0.782983254107, 0. , -0.782983254107 },
|
||||
{ 0. , -0.782983254107, -0.782983254107 },
|
||||
{ 0. , -0.782983254107, 0.782983254107 },
|
||||
{ 0.782983254107, -0.782983254107, 0. },
|
||||
{ 0.782983254107, 0. , 0.782983254107 },
|
||||
{ 0. , 0.782983254107, -0.782983254107 },
|
||||
{ 0.782983254107, 0. , -0.782983254107 },
|
||||
{ 0.782983254107, 0.782983254107, 0. } };
|
||||
|
||||
const double ptm_template_dhex[PTM_NUM_POINTS_DHEX][3] = { { 0. , 0. , 0. },
|
||||
{ -0.391491627053, -0.391491627053, -0.391491627053 },
|
||||
{ 0.391491627053, -0.391491627053, 0.391491627053 },
|
||||
{ -0.391491627053, 0.391491627053, 0.391491627053 },
|
||||
{ 0.391491627053, 0.391491627053, -0.391491627053 },
|
||||
{ -0.260994418036, -1.043977672142, -0.260994418036 },
|
||||
{ -1.043977672142, -0.260994418036, -0.260994418036 },
|
||||
{ -0.260994418036, -0.260994418036, -1.043977672142 },
|
||||
{ 0.782983254107, 0. , 0.782983254107 },
|
||||
{ 0.782983254107, -0.782983254107, 0. },
|
||||
{ 0. , -0.782983254107, 0.782983254107 },
|
||||
{ 0. , 0.782983254107, 0.782983254107 },
|
||||
{ -0.782983254107, 0.782983254107, 0. },
|
||||
{ -0.782983254107, 0. , 0.782983254107 },
|
||||
{ 0.782983254107, 0.782983254107, 0. },
|
||||
{ 0. , 0.782983254107, -0.782983254107 },
|
||||
{ 0.782983254107, 0. , -0.782983254107 } };
|
||||
#endif
|
||||
|
||||
363
src/PTM/ptm_convex_hull_incremental.cpp
Normal file
363
src/PTM/ptm_convex_hull_incremental.cpp
Normal file
@ -0,0 +1,363 @@
|
||||
#include <cmath>
|
||||
#include <cfloat>
|
||||
#include <string.h>
|
||||
#include <cassert>
|
||||
#include <algorithm>
|
||||
#include "ptm_convex_hull_incremental.h"
|
||||
#include "ptm_constants.h"
|
||||
|
||||
|
||||
#define VISIBLE 1
|
||||
#define INVISIBLE 2
|
||||
#define BOTH 3
|
||||
#define TOLERANCE 1E-8
|
||||
|
||||
static double norm_squared(double* p)
|
||||
{
|
||||
double x = p[0];
|
||||
double y = p[1];
|
||||
double z = p[2];
|
||||
|
||||
return x*x + y*y + z*z;
|
||||
}
|
||||
|
||||
static double dot_product(const double* a, const double* b)
|
||||
{
|
||||
return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
|
||||
}
|
||||
|
||||
static void cross_product(double* a, double* b, double* c)
|
||||
{
|
||||
c[0] = a[1] * b[2] - a[2] * b[1];
|
||||
c[1] = a[2] * b[0] - a[0] * b[2];
|
||||
c[2] = a[0] * b[1] - a[1] * b[0];
|
||||
}
|
||||
|
||||
static void calculate_plane_normal(const double (*points)[3], int a, int b, int c, double* plane_normal)
|
||||
{
|
||||
double u[3] = { points[b][0] - points[a][0],
|
||||
points[b][1] - points[a][1],
|
||||
points[b][2] - points[a][2] };
|
||||
|
||||
double v[3] = { points[c][0] - points[a][0],
|
||||
points[c][1] - points[a][1],
|
||||
points[c][2] - points[a][2] };
|
||||
|
||||
cross_product(u, v, plane_normal);
|
||||
double norm = sqrt(norm_squared(plane_normal));
|
||||
plane_normal[0] /= norm;
|
||||
plane_normal[1] /= norm;
|
||||
plane_normal[2] /= norm;
|
||||
}
|
||||
|
||||
static double point_plane_distance(const double* w, const double* plane_point, const double* plane_cross)
|
||||
{
|
||||
return plane_cross[0] * (plane_point[0] - w[0])
|
||||
+ plane_cross[1] * (plane_point[1] - w[1])
|
||||
+ plane_cross[2] * (plane_point[2] - w[2]);
|
||||
}
|
||||
|
||||
static bool calc_max_extent(int num_points, const double (*points)[3], int* min_index, int* max_index)
|
||||
{
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
double dmin = DBL_MAX, dmax = -DBL_MAX;
|
||||
int imin = 0, imax = 0;
|
||||
|
||||
for (int i = 0;i<num_points;i++)
|
||||
{
|
||||
double d = points[i][j];
|
||||
if (d < dmin)
|
||||
{
|
||||
dmin = d;
|
||||
imin = i;
|
||||
}
|
||||
if (d > dmax)
|
||||
{
|
||||
dmax = d;
|
||||
imax = i;
|
||||
}
|
||||
}
|
||||
|
||||
if (imin == imax)
|
||||
return false; //degenerate point set
|
||||
|
||||
min_index[j] = imin;
|
||||
max_index[j] = imax;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
static bool find_third_point(int num_points, const double (*points)[3], int a, int b, int* p_c)
|
||||
{
|
||||
const double* x1 = points[a];
|
||||
const double* x2 = points[b];
|
||||
|
||||
double x2x1[3] = {x2[0] - x1[0], x2[1] - x1[1], x2[2] - x1[2]};
|
||||
double ns_x2x1 = norm_squared(x2x1);
|
||||
|
||||
int bi = -1;
|
||||
double max_dist = 0.0;
|
||||
for (int i = 0;i<num_points;i++)
|
||||
{
|
||||
if (i == a || i == b)
|
||||
continue;
|
||||
|
||||
const double* x0 = points[i];
|
||||
|
||||
double x1x0[3] = {x1[0] - x0[0], x1[1] - x0[1], x1[2] - x0[2]};
|
||||
double dot = dot_product(x1x0, x2x1);
|
||||
double dist = (norm_squared(x1x0) * ns_x2x1 - dot*dot) / ns_x2x1;
|
||||
|
||||
if (dist > max_dist)
|
||||
{
|
||||
max_dist = dist;
|
||||
bi = i;
|
||||
}
|
||||
}
|
||||
|
||||
*p_c = bi;
|
||||
return max_dist > TOLERANCE;
|
||||
}
|
||||
|
||||
static bool find_fourth_point(int num_points, const double (*points)[3], int a, int b, int c, int* p_d)
|
||||
{
|
||||
double plane_normal[3];
|
||||
calculate_plane_normal(points, a, b, c, plane_normal);
|
||||
|
||||
|
||||
int bi = -1;
|
||||
double max_dist = 0.0;
|
||||
for (int i = 0;i<num_points;i++)
|
||||
{
|
||||
if (i == a || i == b || i == c)
|
||||
continue;
|
||||
|
||||
const double* x0 = points[i];
|
||||
double dist = fabs(point_plane_distance(x0, points[a], plane_normal));
|
||||
if (dist > max_dist)
|
||||
{
|
||||
max_dist = dist;
|
||||
bi = i;
|
||||
}
|
||||
}
|
||||
|
||||
*p_d = bi;
|
||||
return max_dist > TOLERANCE;
|
||||
}
|
||||
|
||||
static int initial_simplex(int num_points, const double (*points)[3], int* initial_vertices)
|
||||
{
|
||||
int min_index[3] = {0};
|
||||
int max_index[3] = {0};
|
||||
if (!calc_max_extent(num_points, points, min_index, max_index))
|
||||
return -1;
|
||||
|
||||
int bi = -1;
|
||||
double max_dist = 0.0;
|
||||
for (int i = 0;i<3;i++)
|
||||
{
|
||||
int a = min_index[i], b = max_index[i];
|
||||
double delta[3] = { points[a][0] - points[b][0],
|
||||
points[a][1] - points[b][1],
|
||||
points[a][2] - points[b][2] };
|
||||
double dist = norm_squared(delta);
|
||||
if (dist > max_dist)
|
||||
{
|
||||
bi = i;
|
||||
max_dist = dist;
|
||||
}
|
||||
}
|
||||
|
||||
//first two points are (a, b)
|
||||
int a = min_index[bi], b = max_index[bi], c = -1, d = -1;
|
||||
|
||||
if (!find_third_point(num_points, points, a, b, &c))
|
||||
return -2;
|
||||
|
||||
if (!find_fourth_point(num_points, points, a, b, c, &d))
|
||||
return -3;
|
||||
|
||||
initial_vertices[0] = a;
|
||||
initial_vertices[1] = b;
|
||||
initial_vertices[2] = c;
|
||||
initial_vertices[3] = d;
|
||||
return 0;
|
||||
}
|
||||
|
||||
static bool visible(const double* w, const double* plane_point, const double* plane_normal)
|
||||
{
|
||||
return point_plane_distance(w, plane_point, plane_normal) > 0;
|
||||
}
|
||||
|
||||
void add_facet(const double (*points)[3], int a, int b, int c, int8_t* facet, double* plane_normal, double* barycentre)
|
||||
{
|
||||
calculate_plane_normal(points, a, b, c, plane_normal);
|
||||
if (visible(barycentre, points[a], plane_normal))
|
||||
{
|
||||
plane_normal[0] = -plane_normal[0];
|
||||
plane_normal[1] = -plane_normal[1];
|
||||
plane_normal[2] = -plane_normal[2];
|
||||
|
||||
facet[0] = b;
|
||||
facet[1] = a;
|
||||
facet[2] = c;
|
||||
}
|
||||
else
|
||||
{
|
||||
facet[0] = a;
|
||||
facet[1] = b;
|
||||
facet[2] = c;
|
||||
}
|
||||
}
|
||||
|
||||
static int initialize_convex_hull(int num_points, const double (*points)[3], int8_t facets[][3], double plane_normal[][3], bool* processed, int* initial_vertices, double* barycentre)
|
||||
{
|
||||
memset(processed, 0, PTM_MAX_POINTS * sizeof(bool));
|
||||
memset(barycentre, 0, 3 * sizeof(double));
|
||||
int ret = initial_simplex(num_points, points, initial_vertices);
|
||||
if (ret != 0)
|
||||
return ret;
|
||||
|
||||
for (int i = 0;i<4;i++)
|
||||
{
|
||||
int a = initial_vertices[i];
|
||||
processed[a] = true;
|
||||
|
||||
barycentre[0] += points[a][0];
|
||||
barycentre[1] += points[a][1];
|
||||
barycentre[2] += points[a][2];
|
||||
}
|
||||
barycentre[0] /= 4;
|
||||
barycentre[1] /= 4;
|
||||
barycentre[2] /= 4;
|
||||
|
||||
add_facet(points, initial_vertices[0], initial_vertices[1], initial_vertices[2], facets[0], plane_normal[0], barycentre);
|
||||
add_facet(points, initial_vertices[0], initial_vertices[1], initial_vertices[3], facets[1], plane_normal[1], barycentre);
|
||||
add_facet(points, initial_vertices[0], initial_vertices[2], initial_vertices[3], facets[2], plane_normal[2], barycentre);
|
||||
add_facet(points, initial_vertices[1], initial_vertices[2], initial_vertices[3], facets[3], plane_normal[3], barycentre);
|
||||
return 0;
|
||||
}
|
||||
|
||||
int get_convex_hull(int num_points, const double (*points)[3], convexhull_t* ch, int8_t simplex[][3])
|
||||
{
|
||||
assert( num_points == PTM_NUM_POINTS_FCC
|
||||
|| num_points == PTM_NUM_POINTS_HCP
|
||||
|| num_points == PTM_NUM_POINTS_BCC
|
||||
|| num_points == PTM_NUM_POINTS_ICO
|
||||
|| num_points == PTM_NUM_POINTS_SC
|
||||
|| num_points == PTM_NUM_POINTS_DCUB
|
||||
|| num_points == PTM_NUM_POINTS_DHEX);
|
||||
|
||||
int ret = 0;
|
||||
int num_prev = ch->num_prev;
|
||||
ch->num_prev = num_points;
|
||||
if (!ch->ok || 0)
|
||||
{
|
||||
ret = initialize_convex_hull(num_points, points, ch->facets, ch->plane_normal, ch->processed, ch->initial_vertices, ch->barycentre);
|
||||
if (ret != 0)
|
||||
return ret;
|
||||
|
||||
ch->num_facets = 4;
|
||||
num_prev = 0;
|
||||
}
|
||||
|
||||
for (int i = num_prev;i<num_points;i++)
|
||||
{
|
||||
if (ch->processed[i])
|
||||
continue;
|
||||
ch->processed[i] = true;
|
||||
|
||||
int num_to_add = 0;
|
||||
int8_t to_add[PTM_MAX_FACETS][3];
|
||||
int8_t edge_visible[PTM_MAX_POINTS][PTM_MAX_POINTS];
|
||||
memset(edge_visible, 0, sizeof(int8_t) * PTM_MAX_POINTS * PTM_MAX_POINTS);
|
||||
for (int j = 0;j<ch->num_facets;j++)
|
||||
{
|
||||
int a = ch->facets[j][0];
|
||||
int b = ch->facets[j][1];
|
||||
int c = ch->facets[j][2];
|
||||
|
||||
int u = 0, v = 0, w = 0;
|
||||
|
||||
double distance = point_plane_distance(points[i], points[a], ch->plane_normal[j]);
|
||||
bool vis = distance > TOLERANCE;
|
||||
if (vis)
|
||||
{
|
||||
u = edge_visible[a][b] |= VISIBLE;
|
||||
edge_visible[b][a] |= VISIBLE;
|
||||
|
||||
v = edge_visible[b][c] |= VISIBLE;
|
||||
edge_visible[c][b] |= VISIBLE;
|
||||
|
||||
w = edge_visible[c][a] |= VISIBLE;
|
||||
edge_visible[a][c] |= VISIBLE;
|
||||
|
||||
memcpy(ch->facets[j], ch->facets[ch->num_facets-1], 3 * sizeof(int8_t));
|
||||
memcpy(ch->plane_normal[j], ch->plane_normal[ch->num_facets-1], 3 * sizeof(double));
|
||||
ch->num_facets--;
|
||||
j--;
|
||||
}
|
||||
else
|
||||
{
|
||||
u = edge_visible[a][b] |= INVISIBLE;
|
||||
edge_visible[b][a] |= INVISIBLE;
|
||||
|
||||
v = edge_visible[b][c] |= INVISIBLE;
|
||||
edge_visible[c][b] |= INVISIBLE;
|
||||
|
||||
w = edge_visible[c][a] |= INVISIBLE;
|
||||
edge_visible[a][c] |= INVISIBLE;
|
||||
}
|
||||
|
||||
if (u == BOTH)
|
||||
{
|
||||
to_add[num_to_add][0] = i;
|
||||
to_add[num_to_add][1] = a;
|
||||
to_add[num_to_add][2] = b;
|
||||
num_to_add++;
|
||||
}
|
||||
|
||||
if (v == BOTH)
|
||||
{
|
||||
to_add[num_to_add][0] = i;
|
||||
to_add[num_to_add][1] = b;
|
||||
to_add[num_to_add][2] = c;
|
||||
num_to_add++;
|
||||
}
|
||||
|
||||
if (w == BOTH)
|
||||
{
|
||||
to_add[num_to_add][0] = i;
|
||||
to_add[num_to_add][1] = c;
|
||||
to_add[num_to_add][2] = a;
|
||||
num_to_add++;
|
||||
}
|
||||
}
|
||||
|
||||
for (int j = 0;j<num_to_add;j++)
|
||||
{
|
||||
if (ch->num_facets >= PTM_MAX_FACETS)
|
||||
return -4;
|
||||
|
||||
add_facet(points, to_add[j][0], to_add[j][1], to_add[j][2], ch->facets[ch->num_facets], ch->plane_normal[ch->num_facets], ch->barycentre); ch->num_facets++;
|
||||
}
|
||||
}
|
||||
|
||||
for (int i=0;i<ch->num_facets;i++)
|
||||
{
|
||||
int a = ch->facets[i][0];
|
||||
int b = ch->facets[i][1];
|
||||
int c = ch->facets[i][2];
|
||||
if (a == 0 || b == 0 || c == 0)
|
||||
return 1; //central atom contained in convex hull
|
||||
|
||||
simplex[i][0] = a - 1;
|
||||
simplex[i][1] = b - 1;
|
||||
simplex[i][2] = c - 1;
|
||||
}
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
27
src/PTM/ptm_convex_hull_incremental.h
Normal file
27
src/PTM/ptm_convex_hull_incremental.h
Normal file
@ -0,0 +1,27 @@
|
||||
#ifndef PTM_CONVEX_HULL_INCREMENTAL_H
|
||||
#define PTM_CONVEX_HULL_INCREMENTAL_H
|
||||
|
||||
|
||||
#include <stdint.h>
|
||||
#include <stdbool.h>
|
||||
#include "ptm_constants.h"
|
||||
|
||||
|
||||
typedef struct
|
||||
{
|
||||
int8_t facets[PTM_MAX_FACETS][3];
|
||||
double plane_normal[PTM_MAX_FACETS][3];
|
||||
bool processed[PTM_MAX_POINTS];
|
||||
int initial_vertices[4];
|
||||
double barycentre[3];
|
||||
int num_facets;
|
||||
int num_prev;
|
||||
bool ok;
|
||||
|
||||
} convexhull_t;
|
||||
|
||||
void add_facet(const double (*points)[3], int a, int b, int c, int8_t* facet, double* plane_normal, double* barycentre);
|
||||
int get_convex_hull(int num_points, const double (*points)[3], convexhull_t* ch, int8_t simplex[][3]);
|
||||
|
||||
#endif
|
||||
|
||||
37
src/PTM/ptm_deformation_gradient.cpp
Normal file
37
src/PTM/ptm_deformation_gradient.cpp
Normal file
@ -0,0 +1,37 @@
|
||||
#include "ptm_deformation_gradient.h"
|
||||
|
||||
|
||||
void calculate_deformation_gradient(int num_points, const double (*ideal_points)[3], int8_t* mapping, double (*normalized)[3], const double (*penrose)[3], double* F, double* res)
|
||||
{
|
||||
for (int i = 0;i<3;i++)
|
||||
{
|
||||
for (int j = 0;j<3;j++)
|
||||
{
|
||||
double acc = 0.0;
|
||||
for (int k = 0;k<num_points;k++)
|
||||
acc += penrose[k][j] * normalized[mapping[k]][i];
|
||||
|
||||
F[i*3 + j] = acc;
|
||||
}
|
||||
}
|
||||
|
||||
res[0] = 0;
|
||||
res[1] = 0;
|
||||
res[2] = 0;
|
||||
|
||||
for (int k = 0;k<num_points;k++)
|
||||
{
|
||||
for (int i = 0;i<3;i++)
|
||||
{
|
||||
double acc = 0.0;
|
||||
for (int j = 0;j<3;j++)
|
||||
{
|
||||
acc += F[i*3 + j] * ideal_points[k][j];
|
||||
}
|
||||
|
||||
double delta = acc - normalized[mapping[k]][i];
|
||||
res[i] += delta * delta;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
142
src/PTM/ptm_deformation_gradient.h
Normal file
142
src/PTM/ptm_deformation_gradient.h
Normal file
@ -0,0 +1,142 @@
|
||||
#ifndef PTM_DEFORMATION_GRADIENT_H
|
||||
#define PTM_DEFORMATION_GRADIENT_H
|
||||
|
||||
#include <stdint.h>
|
||||
#include "ptm_constants.h"
|
||||
|
||||
void calculate_deformation_gradient(int num_points, const double (*ideal_points)[3], int8_t* mapping, double (*normalized)[3], const double (*penrose)[3], double* F, double* res);
|
||||
|
||||
//sc
|
||||
#define k_sc 0.5
|
||||
const double penrose_sc[PTM_NUM_POINTS_SC][3] = {
|
||||
{0, 0, 0},
|
||||
{0, 0, -k_sc},
|
||||
{0, 0, k_sc},
|
||||
{0, -k_sc, 0},
|
||||
{0, k_sc, 0},
|
||||
{-k_sc, 0, 0},
|
||||
{k_sc, 0, 0},
|
||||
};
|
||||
|
||||
//fcc
|
||||
#define k_fcc 0.17677669529663678216
|
||||
const double penrose_fcc[PTM_NUM_POINTS_FCC][3] = {
|
||||
{0, 0, 0},
|
||||
{0, k_fcc, k_fcc},
|
||||
{0, -k_fcc, -k_fcc},
|
||||
{0, k_fcc, -k_fcc},
|
||||
{0, -k_fcc, k_fcc},
|
||||
{k_fcc, 0, k_fcc},
|
||||
{-k_fcc, 0, -k_fcc},
|
||||
{k_fcc, 0, -k_fcc},
|
||||
{-k_fcc, 0, k_fcc},
|
||||
{k_fcc, k_fcc, -0},
|
||||
{-k_fcc, -k_fcc, 0},
|
||||
{k_fcc, -k_fcc, 0},
|
||||
{-k_fcc, k_fcc, -0},
|
||||
};
|
||||
|
||||
//hcp
|
||||
#define k_hcp 0.17677669529663678216
|
||||
const double penrose_hcp[PTM_NUM_POINTS_HCP][3] = {
|
||||
{0, 0, 0},
|
||||
{k_hcp, 0, k_hcp},
|
||||
{-k_hcp/3, -4*k_hcp/3, -k_hcp/3},
|
||||
{k_hcp, k_hcp, 0},
|
||||
{-k_hcp/3, -k_hcp/3, -4*k_hcp/3},
|
||||
{0, k_hcp, k_hcp},
|
||||
{-4*k_hcp/3, -k_hcp/3, -k_hcp/3},
|
||||
{-k_hcp, k_hcp, -0},
|
||||
{0, k_hcp, -k_hcp},
|
||||
{k_hcp, 0, -k_hcp},
|
||||
{k_hcp, -k_hcp, 0},
|
||||
{-k_hcp, 0, k_hcp},
|
||||
{0, -k_hcp, k_hcp},
|
||||
};
|
||||
|
||||
//ico
|
||||
#define k_ico 0.13143277802974323576
|
||||
#define phi 1.61803398874989490253
|
||||
//((1.0 + sqrt(5)) / 2)
|
||||
const double penrose_ico[PTM_NUM_POINTS_ICO][3] = {
|
||||
{0, 0, 0},
|
||||
{0, k_ico, phi*k_ico},
|
||||
{0, -k_ico, -phi*k_ico},
|
||||
{0, k_ico, -phi*k_ico},
|
||||
{0, -k_ico, phi*k_ico},
|
||||
{-k_ico, -phi*k_ico, -0},
|
||||
{k_ico, phi*k_ico, 0},
|
||||
{k_ico, -phi*k_ico, 0},
|
||||
{-k_ico, phi*k_ico, -0},
|
||||
{-phi*k_ico, 0, -k_ico},
|
||||
{phi*k_ico, 0, k_ico},
|
||||
{phi*k_ico, 0, -k_ico},
|
||||
{-phi*k_ico, 0, k_ico},
|
||||
};
|
||||
|
||||
//bcc
|
||||
#define k_bcc 0.11543038598460284017
|
||||
const double penrose_bcc[PTM_NUM_POINTS_BCC][3] = {
|
||||
{0, 0, 0},
|
||||
{-k_bcc, -k_bcc, -k_bcc},
|
||||
{k_bcc, k_bcc, k_bcc},
|
||||
{k_bcc, -k_bcc, -k_bcc},
|
||||
{-k_bcc, k_bcc, k_bcc},
|
||||
{-k_bcc, k_bcc, -k_bcc},
|
||||
{k_bcc, -k_bcc, k_bcc},
|
||||
{-k_bcc, -k_bcc, k_bcc},
|
||||
{k_bcc, k_bcc, -k_bcc},
|
||||
{0, 0, -2*k_bcc},
|
||||
{0, 0, 2*k_bcc},
|
||||
{0, -2*k_bcc, 0},
|
||||
{0, 2*k_bcc, 0},
|
||||
{-2*k_bcc, 0, 0},
|
||||
{2*k_bcc, 0, -0},
|
||||
};
|
||||
|
||||
//dcub
|
||||
#define kdcub 0.07095369570691034689
|
||||
const double penrose_dcub[PTM_NUM_POINTS_DCUB][3] = {
|
||||
{ 0, 0, 0 },
|
||||
{ -kdcub, kdcub, kdcub },
|
||||
{ -kdcub, -kdcub, -kdcub },
|
||||
{ kdcub, -kdcub, kdcub },
|
||||
{ kdcub, kdcub, -kdcub },
|
||||
{ -2 * kdcub, 0, 2 * kdcub },
|
||||
{ -2 * kdcub, 2 * kdcub, 0 },
|
||||
{ 0, 2 * kdcub, 2 * kdcub },
|
||||
{ -2 * kdcub, -2 * kdcub, 0 },
|
||||
{ -2 * kdcub, 0, -2 * kdcub },
|
||||
{ 0, -2 * kdcub, -2 * kdcub },
|
||||
{ 0, -2 * kdcub, 2 * kdcub },
|
||||
{ 2 * kdcub, -2 * kdcub, 0 },
|
||||
{ 2 * kdcub, 0, 2 * kdcub },
|
||||
{ 0, 2 * kdcub, -2 * kdcub },
|
||||
{ 2 * kdcub, 0, -2 * kdcub },
|
||||
{ 2 * kdcub, 2 * kdcub, 0 },
|
||||
};
|
||||
|
||||
|
||||
#define kdhex 0.04730246380471011397
|
||||
const double penrose_dhex[PTM_NUM_POINTS_DHEX][3] = {
|
||||
{ 0, 0, 0 },
|
||||
{ -kdcub, -kdcub, -kdcub },
|
||||
{ kdcub, -kdcub, kdcub },
|
||||
{ -kdcub, kdcub, kdcub },
|
||||
{ kdcub, kdcub, -kdcub },
|
||||
{ -kdhex, -4 * kdhex, -kdhex },
|
||||
{ -4 * kdhex, -kdhex, -kdhex },
|
||||
{ -kdhex, -kdhex, -4 * kdhex },
|
||||
{ 2 * kdcub, 0, 2 * kdcub },
|
||||
{ 2 * kdcub, -2 * kdcub, 0 },
|
||||
{ 0, -2 * kdcub, 2 * kdcub },
|
||||
{ 0, 2 * kdcub, 2 * kdcub },
|
||||
{ -2 * kdcub, 2 * kdcub, 0 },
|
||||
{ -2 * kdcub, 0, 2 * kdcub },
|
||||
{ 2 * kdcub, 2 * kdcub, 0 },
|
||||
{ 0, 2 * kdcub, -2 * kdcub },
|
||||
{ 2 * kdcub, 0, -2 * kdcub },
|
||||
};
|
||||
#endif
|
||||
|
||||
|
||||
27
src/PTM/ptm_functions.h
Normal file
27
src/PTM/ptm_functions.h
Normal file
@ -0,0 +1,27 @@
|
||||
#ifndef PTM_FUNCTIONS_H
|
||||
#define PTM_FUNCTIONS_H
|
||||
|
||||
#include <stdint.h>
|
||||
#include <stdbool.h>
|
||||
#include "ptm_initialize_data.h"
|
||||
#include "ptm_constants.h"
|
||||
|
||||
|
||||
//------------------------------------
|
||||
// function declarations
|
||||
//------------------------------------
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
|
||||
int ptm_index( ptm_local_handle_t local_handle, int32_t flags, int num_points, double (*atomic_positions)[3], int32_t* atomic_numbers, bool topological_ordering, //inputs
|
||||
int32_t* p_type, int32_t* p_alloy_type, double* p_scale, double* p_rmsd, double* q, double* F, double* F_res, double* U, double* P, int8_t* mapping, double* p_interatomic_distance, double* p_lattice_constant); //outputs
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
||||
180
src/PTM/ptm_fundamental_mappings.h
Normal file
180
src/PTM/ptm_fundamental_mappings.h
Normal file
@ -0,0 +1,180 @@
|
||||
#ifndef PTM_FUNDAMENTAL_MAPPINGS_H
|
||||
#define PTM_FUNDAMENTAL_MAPPINGS_H
|
||||
|
||||
#include <stdint.h>
|
||||
|
||||
#define NUM_CUBIC_MAPPINGS 24
|
||||
#define NUM_ICO_MAPPINGS 60
|
||||
#define NUM_HEX_MAPPINGS 6
|
||||
#define NUM_DCUB_MAPPINGS 12
|
||||
#define NUM_DHEX_MAPPINGS 3
|
||||
|
||||
const int8_t mapping_sc[NUM_CUBIC_MAPPINGS][PTM_MAX_POINTS] = {
|
||||
{0, 1, 2, 3, 4, 5, 6},
|
||||
{0, 2, 1, 4, 3, 5, 6},
|
||||
{0, 2, 1, 3, 4, 6, 5},
|
||||
{0, 1, 2, 4, 3, 6, 5},
|
||||
{0, 3, 4, 5, 6, 1, 2},
|
||||
{0, 5, 6, 2, 1, 4, 3},
|
||||
{0, 6, 5, 1, 2, 4, 3},
|
||||
{0, 4, 3, 5, 6, 2, 1},
|
||||
{0, 5, 6, 1, 2, 3, 4},
|
||||
{0, 4, 3, 6, 5, 1, 2},
|
||||
{0, 3, 4, 6, 5, 2, 1},
|
||||
{0, 6, 5, 2, 1, 3, 4},
|
||||
{0, 3, 4, 2, 1, 5, 6},
|
||||
{0, 6, 5, 3, 4, 1, 2},
|
||||
{0, 1, 2, 5, 6, 4, 3},
|
||||
{0, 4, 3, 1, 2, 5, 6},
|
||||
{0, 5, 6, 3, 4, 2, 1},
|
||||
{0, 1, 2, 6, 5, 3, 4},
|
||||
{0, 2, 1, 5, 6, 3, 4},
|
||||
{0, 5, 6, 4, 3, 1, 2},
|
||||
{0, 3, 4, 1, 2, 6, 5},
|
||||
{0, 2, 1, 6, 5, 4, 3},
|
||||
{0, 6, 5, 4, 3, 2, 1},
|
||||
{0, 4, 3, 2, 1, 6, 5} };
|
||||
|
||||
const int8_t mapping_fcc[NUM_CUBIC_MAPPINGS][PTM_MAX_POINTS] = {
|
||||
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12},
|
||||
{0, 2, 1, 4, 3, 7, 8, 5, 6, 11, 12, 9, 10},
|
||||
{0, 3, 4, 1, 2, 6, 5, 8, 7, 12, 11, 10, 9},
|
||||
{0, 4, 3, 2, 1, 8, 7, 6, 5, 10, 9, 12, 11},
|
||||
{0, 9, 10, 11, 12, 1, 2, 4, 3, 5, 6, 8, 7},
|
||||
{0, 7, 8, 6, 5, 11, 12, 10, 9, 2, 1, 4, 3},
|
||||
{0, 8, 7, 5, 6, 10, 9, 11, 12, 4, 3, 2, 1},
|
||||
{0, 11, 12, 9, 10, 2, 1, 3, 4, 7, 8, 6, 5},
|
||||
{0, 5, 6, 8, 7, 9, 10, 12, 11, 1, 2, 3, 4},
|
||||
{0, 10, 9, 12, 11, 4, 3, 1, 2, 8, 7, 5, 6},
|
||||
{0, 12, 11, 10, 9, 3, 4, 2, 1, 6, 5, 7, 8},
|
||||
{0, 6, 5, 7, 8, 12, 11, 9, 10, 3, 4, 1, 2},
|
||||
{0, 3, 4, 2, 1, 9, 10, 11, 12, 7, 8, 5, 6},
|
||||
{0, 12, 11, 9, 10, 8, 7, 5, 6, 1, 2, 4, 3},
|
||||
{0, 5, 6, 7, 8, 4, 3, 2, 1, 11, 12, 10, 9},
|
||||
{0, 4, 3, 1, 2, 11, 12, 9, 10, 5, 6, 7, 8},
|
||||
{0, 9, 10, 12, 11, 7, 8, 6, 5, 3, 4, 2, 1},
|
||||
{0, 8, 7, 6, 5, 1, 2, 3, 4, 12, 11, 9, 10},
|
||||
{0, 7, 8, 5, 6, 3, 4, 1, 2, 9, 10, 12, 11},
|
||||
{0, 11, 12, 10, 9, 5, 6, 8, 7, 4, 3, 1, 2},
|
||||
{0, 1, 2, 4, 3, 12, 11, 10, 9, 8, 7, 6, 5},
|
||||
{0, 6, 5, 8, 7, 2, 1, 4, 3, 10, 9, 11, 12},
|
||||
{0, 10, 9, 11, 12, 6, 5, 7, 8, 2, 1, 3, 4},
|
||||
{0, 2, 1, 3, 4, 10, 9, 12, 11, 6, 5, 8, 7} };
|
||||
|
||||
const int8_t mapping_bcc[NUM_CUBIC_MAPPINGS][PTM_MAX_POINTS] = {
|
||||
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14},
|
||||
{0, 4, 3, 2, 1, 7, 8, 5, 6, 10, 9, 12, 11, 13, 14},
|
||||
{0, 6, 5, 7, 8, 2, 1, 3, 4, 10, 9, 11, 12, 14, 13},
|
||||
{0, 8, 7, 5, 6, 3, 4, 2, 1, 9, 10, 12, 11, 14, 13},
|
||||
{0, 1, 2, 7, 8, 3, 4, 5, 6, 11, 12, 13, 14, 9, 10},
|
||||
{0, 4, 3, 7, 8, 5, 6, 2, 1, 13, 14, 10, 9, 12, 11},
|
||||
{0, 8, 7, 3, 4, 2, 1, 5, 6, 14, 13, 9, 10, 12, 11},
|
||||
{0, 4, 3, 5, 6, 2, 1, 7, 8, 12, 11, 13, 14, 10, 9},
|
||||
{0, 1, 2, 5, 6, 7, 8, 3, 4, 13, 14, 9, 10, 11, 12},
|
||||
{0, 8, 7, 2, 1, 5, 6, 3, 4, 12, 11, 14, 13, 9, 10},
|
||||
{0, 6, 5, 3, 4, 7, 8, 2, 1, 11, 12, 14, 13, 10, 9},
|
||||
{0, 6, 5, 2, 1, 3, 4, 7, 8, 14, 13, 10, 9, 11, 12},
|
||||
{0, 7, 8, 6, 5, 1, 2, 4, 3, 11, 12, 10, 9, 13, 14},
|
||||
{0, 3, 4, 6, 5, 8, 7, 1, 2, 14, 13, 11, 12, 9, 10},
|
||||
{0, 5, 6, 1, 2, 8, 7, 4, 3, 9, 10, 13, 14, 12, 11},
|
||||
{0, 5, 6, 8, 7, 4, 3, 1, 2, 12, 11, 9, 10, 13, 14},
|
||||
{0, 7, 8, 1, 2, 4, 3, 6, 5, 13, 14, 11, 12, 10, 9},
|
||||
{0, 3, 4, 8, 7, 1, 2, 6, 5, 9, 10, 14, 13, 11, 12},
|
||||
{0, 7, 8, 4, 3, 6, 5, 1, 2, 10, 9, 13, 14, 11, 12},
|
||||
{0, 5, 6, 4, 3, 1, 2, 8, 7, 13, 14, 12, 11, 9, 10},
|
||||
{0, 3, 4, 1, 2, 6, 5, 8, 7, 11, 12, 9, 10, 14, 13},
|
||||
{0, 2, 1, 6, 5, 4, 3, 8, 7, 10, 9, 14, 13, 12, 11},
|
||||
{0, 2, 1, 8, 7, 6, 5, 4, 3, 14, 13, 12, 11, 10, 9},
|
||||
{0, 2, 1, 4, 3, 8, 7, 6, 5, 12, 11, 10, 9, 14, 13} };
|
||||
|
||||
const int8_t mapping_ico[NUM_ICO_MAPPINGS][PTM_MAX_POINTS] = {
|
||||
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12},
|
||||
{0, 10, 9, 8, 7, 5, 6, 2, 1, 12, 11, 3, 4},
|
||||
{0, 1, 2, 9, 10, 7, 8, 11, 12, 5, 6, 3, 4},
|
||||
{0, 4, 3, 8, 7, 2, 1, 11, 12, 9, 10, 6, 5},
|
||||
{0, 6, 5, 9, 10, 4, 3, 7, 8, 12, 11, 2, 1},
|
||||
{0, 12, 11, 3, 4, 7, 8, 10, 9, 2, 1, 6, 5},
|
||||
{0, 4, 3, 6, 5, 9, 10, 2, 1, 8, 7, 11, 12},
|
||||
{0, 8, 7, 2, 1, 4, 3, 10, 9, 5, 6, 11, 12},
|
||||
{0, 10, 9, 3, 4, 12, 11, 5, 6, 8, 7, 2, 1},
|
||||
{0, 12, 11, 6, 5, 2, 1, 7, 8, 3, 4, 10, 9},
|
||||
{0, 1, 2, 11, 12, 9, 10, 5, 6, 3, 4, 7, 8},
|
||||
{0, 8, 7, 11, 12, 5, 6, 4, 3, 2, 1, 10, 9},
|
||||
{0, 6, 5, 2, 1, 12, 11, 4, 3, 9, 10, 7, 8},
|
||||
{0, 3, 4, 5, 6, 1, 2, 10, 9, 12, 11, 7, 8},
|
||||
{0, 3, 4, 7, 8, 12, 11, 1, 2, 5, 6, 10, 9},
|
||||
{0, 6, 5, 7, 8, 9, 10, 12, 11, 2, 1, 4, 3},
|
||||
{0, 9, 10, 11, 12, 4, 3, 1, 2, 7, 8, 6, 5},
|
||||
{0, 11, 12, 9, 10, 1, 2, 4, 3, 8, 7, 5, 6},
|
||||
{0, 8, 7, 5, 6, 10, 9, 11, 12, 4, 3, 2, 1},
|
||||
{0, 10, 9, 2, 1, 8, 7, 12, 11, 3, 4, 5, 6},
|
||||
{0, 12, 11, 2, 1, 10, 9, 6, 5, 7, 8, 3, 4},
|
||||
{0, 9, 10, 6, 5, 7, 8, 4, 3, 11, 12, 1, 2},
|
||||
{0, 8, 7, 10, 9, 2, 1, 5, 6, 11, 12, 4, 3},
|
||||
{0, 6, 5, 12, 11, 7, 8, 2, 1, 4, 3, 9, 10},
|
||||
{0, 11, 12, 8, 7, 4, 3, 5, 6, 1, 2, 9, 10},
|
||||
{0, 4, 3, 11, 12, 8, 7, 9, 10, 6, 5, 2, 1},
|
||||
{0, 4, 3, 9, 10, 11, 12, 6, 5, 2, 1, 8, 7},
|
||||
{0, 12, 11, 10, 9, 3, 4, 2, 1, 6, 5, 7, 8},
|
||||
{0, 5, 6, 8, 7, 11, 12, 10, 9, 3, 4, 1, 2},
|
||||
{0, 7, 8, 6, 5, 12, 11, 9, 10, 1, 2, 3, 4},
|
||||
{0, 10, 9, 12, 11, 2, 1, 3, 4, 5, 6, 8, 7},
|
||||
{0, 7, 8, 1, 2, 9, 10, 3, 4, 12, 11, 6, 5},
|
||||
{0, 5, 6, 1, 2, 3, 4, 11, 12, 8, 7, 10, 9},
|
||||
{0, 7, 8, 12, 11, 3, 4, 6, 5, 9, 10, 1, 2},
|
||||
{0, 1, 2, 5, 6, 11, 12, 3, 4, 7, 8, 9, 10},
|
||||
{0, 11, 12, 1, 2, 5, 6, 9, 10, 4, 3, 8, 7},
|
||||
{0, 5, 6, 3, 4, 10, 9, 1, 2, 11, 12, 8, 7},
|
||||
{0, 5, 6, 10, 9, 8, 7, 3, 4, 1, 2, 11, 12},
|
||||
{0, 3, 4, 12, 11, 10, 9, 7, 8, 1, 2, 5, 6},
|
||||
{0, 9, 10, 7, 8, 1, 2, 6, 5, 4, 3, 11, 12},
|
||||
{0, 9, 10, 1, 2, 11, 12, 7, 8, 6, 5, 4, 3},
|
||||
{0, 7, 8, 3, 4, 1, 2, 12, 11, 6, 5, 9, 10},
|
||||
{0, 11, 12, 5, 6, 8, 7, 1, 2, 9, 10, 4, 3},
|
||||
{0, 1, 2, 7, 8, 3, 4, 9, 10, 11, 12, 5, 6},
|
||||
{0, 3, 4, 10, 9, 5, 6, 12, 11, 7, 8, 1, 2},
|
||||
{0, 2, 1, 4, 3, 8, 7, 6, 5, 12, 11, 10, 9},
|
||||
{0, 2, 1, 12, 11, 6, 5, 10, 9, 8, 7, 4, 3},
|
||||
{0, 9, 10, 4, 3, 6, 5, 11, 12, 1, 2, 7, 8},
|
||||
{0, 11, 12, 4, 3, 9, 10, 8, 7, 5, 6, 1, 2},
|
||||
{0, 2, 1, 10, 9, 12, 11, 8, 7, 4, 3, 6, 5},
|
||||
{0, 5, 6, 11, 12, 1, 2, 8, 7, 10, 9, 3, 4},
|
||||
{0, 10, 9, 5, 6, 3, 4, 8, 7, 2, 1, 12, 11},
|
||||
{0, 12, 11, 7, 8, 6, 5, 3, 4, 10, 9, 2, 1},
|
||||
{0, 7, 8, 9, 10, 6, 5, 1, 2, 3, 4, 12, 11},
|
||||
{0, 2, 1, 8, 7, 10, 9, 4, 3, 6, 5, 12, 11},
|
||||
{0, 8, 7, 4, 3, 11, 12, 2, 1, 10, 9, 5, 6},
|
||||
{0, 6, 5, 4, 3, 2, 1, 9, 10, 7, 8, 12, 11},
|
||||
{0, 2, 1, 6, 5, 4, 3, 12, 11, 10, 9, 8, 7},
|
||||
{0, 3, 4, 1, 2, 7, 8, 5, 6, 10, 9, 12, 11},
|
||||
{0, 4, 3, 2, 1, 6, 5, 8, 7, 11, 12, 9, 10} };
|
||||
|
||||
const int8_t mapping_hcp[NUM_HEX_MAPPINGS][PTM_MAX_POINTS] = {
|
||||
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12},
|
||||
{0, 5, 6, 1, 2, 3, 4, 9, 10, 12, 11, 8, 7},
|
||||
{0, 3, 4, 5, 6, 1, 2, 12, 11, 7, 8, 10, 9},
|
||||
{0, 4, 3, 2, 1, 6, 5, 11, 12, 10, 9, 7, 8},
|
||||
{0, 2, 1, 6, 5, 4, 3, 8, 7, 11, 12, 9, 10},
|
||||
{0, 6, 5, 4, 3, 2, 1, 10, 9, 8, 7, 12, 11} };
|
||||
|
||||
const int8_t mapping_dcub[NUM_DCUB_MAPPINGS][PTM_MAX_POINTS] = {
|
||||
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16},
|
||||
{0, 2, 1, 4, 3, 9, 8, 10, 6, 5, 7, 14, 16, 15, 11, 13, 12},
|
||||
{0, 4, 3, 2, 1, 15, 16, 14, 12, 13, 11, 10, 8, 9, 7, 5, 6},
|
||||
{0, 3, 4, 1, 2, 13, 12, 11, 16, 15, 14, 7, 6, 5, 10, 9, 8},
|
||||
{0, 4, 2, 1, 3, 14, 15, 16, 9, 10, 8, 6, 5, 7, 12, 11, 13},
|
||||
{0, 4, 1, 3, 2, 16, 14, 15, 7, 6, 5, 13, 11, 12, 9, 8, 10},
|
||||
{0, 1, 4, 2, 3, 6, 7, 5, 14, 16, 15, 9, 10, 8, 13, 12, 11},
|
||||
{0, 3, 1, 2, 4, 11, 13, 12, 5, 7, 6, 8, 9, 10, 16, 14, 15},
|
||||
{0, 3, 2, 4, 1, 12, 11, 13, 10, 8, 9, 15, 14, 16, 5, 6, 7},
|
||||
{0, 2, 4, 3, 1, 10, 9, 8, 15, 14, 16, 12, 13, 11, 6, 7, 5},
|
||||
{0, 1, 3, 4, 2, 7, 5, 6, 13, 11, 12, 16, 15, 14, 8, 10, 9},
|
||||
{0, 2, 3, 1, 4, 8, 10, 9, 11, 12, 13, 5, 7, 6, 15, 16, 14} };
|
||||
|
||||
const int8_t mapping_dhex[NUM_DHEX_MAPPINGS][PTM_MAX_POINTS] = {
|
||||
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16},
|
||||
{0, 1, 3, 4, 2, 6, 7, 5, 11, 13, 12, 14, 16, 15, 8, 9, 10},
|
||||
{0, 1, 4, 2, 3, 7, 5, 6, 14, 15, 16, 8, 10, 9, 11, 13, 12} };
|
||||
|
||||
#endif
|
||||
|
||||
2059
src/PTM/ptm_graph_data.cpp
Normal file
2059
src/PTM/ptm_graph_data.cpp
Normal file
File diff suppressed because it is too large
Load Diff
37
src/PTM/ptm_graph_data.h
Normal file
37
src/PTM/ptm_graph_data.h
Normal file
@ -0,0 +1,37 @@
|
||||
#ifndef PTM_GRAPH_DATA_H
|
||||
#define PTM_GRAPH_DATA_H
|
||||
|
||||
#include <stdint.h>
|
||||
#include "ptm_constants.h"
|
||||
|
||||
|
||||
typedef struct
|
||||
{
|
||||
int id;
|
||||
uint64_t hash;
|
||||
int automorphism_index;
|
||||
int num_automorphisms;
|
||||
int8_t canonical_labelling[PTM_MAX_POINTS];
|
||||
int8_t facets[PTM_MAX_FACETS][3];
|
||||
} graph_t;
|
||||
|
||||
#define NUM_SC_GRAPHS 1
|
||||
#define NUM_ICO_GRAPHS 1
|
||||
#define NUM_FCC_GRAPHS 8
|
||||
#define NUM_HCP_GRAPHS 16
|
||||
#define NUM_BCC_GRAPHS 218
|
||||
#define NUM_DCUB_GRAPHS 12
|
||||
#define NUM_DHEX_GRAPHS 24
|
||||
|
||||
extern int8_t automorphisms[][PTM_MAX_POINTS];
|
||||
|
||||
extern graph_t graphs_sc[NUM_SC_GRAPHS];
|
||||
extern graph_t graphs_fcc[NUM_FCC_GRAPHS];
|
||||
extern graph_t graphs_hcp[NUM_HCP_GRAPHS];
|
||||
extern graph_t graphs_ico[NUM_ICO_GRAPHS];
|
||||
extern graph_t graphs_bcc[NUM_BCC_GRAPHS];
|
||||
extern graph_t graphs_dcub[NUM_DCUB_GRAPHS];
|
||||
extern graph_t graphs_dhex[NUM_DHEX_GRAPHS];
|
||||
|
||||
#endif
|
||||
|
||||
52
src/PTM/ptm_graph_tools.cpp
Normal file
52
src/PTM/ptm_graph_tools.cpp
Normal file
@ -0,0 +1,52 @@
|
||||
#include <string.h>
|
||||
#include <algorithm>
|
||||
#include "ptm_graph_tools.h"
|
||||
#include "ptm_constants.h"
|
||||
|
||||
|
||||
bool build_facet_map(int num_facets, int8_t facets[][3], int8_t common[PTM_MAX_NBRS][PTM_MAX_NBRS])
|
||||
{
|
||||
memset(common, -1, sizeof(int8_t) * PTM_MAX_NBRS * PTM_MAX_NBRS);
|
||||
|
||||
for (int i = 0;i<num_facets;i++)
|
||||
{
|
||||
int a = facets[i][0];
|
||||
int b = facets[i][1];
|
||||
int c = facets[i][2];
|
||||
|
||||
//assert(common[a][b] == -1);
|
||||
//assert(common[b][c] == -1);
|
||||
//assert(common[c][a] == -1);
|
||||
if (common[a][b] != -1 || common[b][c] != -1 || common[c][a] != -1)
|
||||
return false;
|
||||
|
||||
common[a][b] = c;
|
||||
common[b][c] = a;
|
||||
common[c][a] = b;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
int graph_degree(int num_facets, int8_t facets[][3], int num_nodes, int8_t* degree)
|
||||
{
|
||||
memset(degree, 0, sizeof(int8_t) * num_nodes);
|
||||
|
||||
for (int i = 0;i<num_facets;i++)
|
||||
{
|
||||
int a = facets[i][0];
|
||||
int b = facets[i][1];
|
||||
int c = facets[i][2];
|
||||
|
||||
degree[a]++;
|
||||
degree[b]++;
|
||||
degree[c]++;
|
||||
}
|
||||
|
||||
int8_t max_degree = 0;
|
||||
for (int i = 0;i<num_nodes;i++)
|
||||
max_degree = std::max(max_degree, degree[i]);
|
||||
|
||||
return max_degree;
|
||||
}
|
||||
|
||||
11
src/PTM/ptm_graph_tools.h
Normal file
11
src/PTM/ptm_graph_tools.h
Normal file
@ -0,0 +1,11 @@
|
||||
#ifndef PTM_GRAPH_TOOLS_H
|
||||
#define PTM_GRAPH_TOOLS_H
|
||||
|
||||
#include <stdint.h>
|
||||
#include "ptm_constants.h"
|
||||
|
||||
bool build_facet_map(int num_facets, int8_t facets[][3], int8_t common[PTM_MAX_NBRS][PTM_MAX_NBRS]);
|
||||
int graph_degree(int num_facets, int8_t facets[][3], int num_nodes, int8_t* degree);
|
||||
|
||||
#endif
|
||||
|
||||
218
src/PTM/ptm_index.cpp
Normal file
218
src/PTM/ptm_index.cpp
Normal file
@ -0,0 +1,218 @@
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <string.h>
|
||||
#include <cmath>
|
||||
#include <cfloat>
|
||||
#include <cassert>
|
||||
#include <algorithm>
|
||||
#include "ptm_convex_hull_incremental.h"
|
||||
#include "ptm_graph_data.h"
|
||||
#include "ptm_deformation_gradient.h"
|
||||
#include "ptm_alloy_types.h"
|
||||
#include "ptm_neighbour_ordering.h"
|
||||
#include "ptm_normalize_vertices.h"
|
||||
#include "ptm_quat.h"
|
||||
#include "ptm_polar.h"
|
||||
#include "ptm_initialize_data.h"
|
||||
#include "ptm_structure_matcher.h"
|
||||
#include "ptm_functions.h"
|
||||
#include "ptm_constants.h"
|
||||
|
||||
|
||||
//todo: verify that c == norm(template[1])
|
||||
static double calculate_interatomic_distance(int type, double scale)
|
||||
{
|
||||
assert(type >= 1 && type <= 7);
|
||||
double c[8] = {0, 1, 1, (7. - 3.5 * sqrt(3)), 1, 1, sqrt(3) * 4. / (6 * sqrt(2) + sqrt(3)), sqrt(3) * 4. / (6 * sqrt(2) + sqrt(3))};
|
||||
return c[type] / scale;
|
||||
}
|
||||
|
||||
static double calculate_lattice_constant(int type, double interatomic_distance)
|
||||
{
|
||||
assert(type >= 1 && type <= 7);
|
||||
double c[8] = {0, 2 / sqrt(2), 2 / sqrt(2), 2. / sqrt(3), 2 / sqrt(2), 1, 4 / sqrt(3), 4 / sqrt(3)};
|
||||
return c[type] * interatomic_distance;
|
||||
}
|
||||
|
||||
static int rotate_into_fundamental_zone(int type, double* q)
|
||||
{
|
||||
if (type == PTM_MATCH_SC) return rotate_quaternion_into_cubic_fundamental_zone(q);
|
||||
if (type == PTM_MATCH_FCC) return rotate_quaternion_into_cubic_fundamental_zone(q);
|
||||
if (type == PTM_MATCH_BCC) return rotate_quaternion_into_cubic_fundamental_zone(q);
|
||||
if (type == PTM_MATCH_ICO) return rotate_quaternion_into_icosahedral_fundamental_zone(q);
|
||||
if (type == PTM_MATCH_HCP) return rotate_quaternion_into_hcp_fundamental_zone(q);
|
||||
if (type == PTM_MATCH_DCUB) return rotate_quaternion_into_diamond_cubic_fundamental_zone(q);
|
||||
if (type == PTM_MATCH_DHEX) return rotate_quaternion_into_diamond_hexagonal_fundamental_zone(q);
|
||||
return -1;
|
||||
}
|
||||
|
||||
static void order_points(ptm_local_handle_t local_handle, int num_points, double (*unpermuted_points)[3], int32_t* unpermuted_numbers, bool topological_ordering,
|
||||
int8_t* ordering, double (*points)[3], int32_t* numbers)
|
||||
{
|
||||
if (topological_ordering)
|
||||
{
|
||||
double normalized_points[PTM_MAX_INPUT_POINTS][3];
|
||||
normalize_vertices(num_points, unpermuted_points, normalized_points);
|
||||
int ret = calculate_neighbour_ordering((void*)local_handle, num_points, (const double (*)[3])normalized_points, ordering);
|
||||
if (ret != 0)
|
||||
topological_ordering = false;
|
||||
}
|
||||
|
||||
if (!topological_ordering)
|
||||
for (int i=0;i<num_points;i++)
|
||||
ordering[i] = i;
|
||||
|
||||
for (int i=0;i<num_points;i++)
|
||||
{
|
||||
memcpy(points[i], &unpermuted_points[ordering[i]], 3 * sizeof(double));
|
||||
|
||||
if (unpermuted_numbers != NULL)
|
||||
numbers[i] = unpermuted_numbers[ordering[i]];
|
||||
}
|
||||
}
|
||||
|
||||
static void output_data(result_t* res, int num_points, int32_t* unpermuted_numbers, double (*points)[3], int32_t* numbers, int8_t* ordering,
|
||||
int32_t* p_type, int32_t* p_alloy_type, double* p_scale, double* p_rmsd, double* q, double* F, double* F_res,
|
||||
double* U, double* P, int8_t* mapping, double* p_interatomic_distance, double* p_lattice_constant)
|
||||
{
|
||||
*p_type = PTM_MATCH_NONE;
|
||||
if (p_alloy_type != NULL)
|
||||
*p_alloy_type = PTM_ALLOY_NONE;
|
||||
|
||||
if (mapping != NULL)
|
||||
memset(mapping, -1, num_points * sizeof(int8_t));
|
||||
|
||||
const refdata_t* ref = res->ref_struct;
|
||||
if (ref == NULL)
|
||||
return;
|
||||
|
||||
*p_type = ref->type;
|
||||
if (p_alloy_type != NULL && unpermuted_numbers != NULL)
|
||||
*p_alloy_type = find_alloy_type(ref, res->mapping, numbers);
|
||||
|
||||
int bi = rotate_into_fundamental_zone(ref->type, res->q);
|
||||
int8_t temp[PTM_MAX_POINTS];
|
||||
for (int i=0;i<ref->num_nbrs+1;i++)
|
||||
temp[ref->mapping[bi][i]] = res->mapping[i];
|
||||
|
||||
memcpy(res->mapping, temp, (ref->num_nbrs+1) * sizeof(int8_t));
|
||||
|
||||
if (F != NULL && F_res != NULL)
|
||||
{
|
||||
double scaled_points[PTM_MAX_INPUT_POINTS][3];
|
||||
|
||||
subtract_barycentre(ref->num_nbrs + 1, points, scaled_points);
|
||||
for (int i = 0;i<ref->num_nbrs + 1;i++)
|
||||
{
|
||||
scaled_points[i][0] *= res->scale;
|
||||
scaled_points[i][1] *= res->scale;
|
||||
scaled_points[i][2] *= res->scale;
|
||||
}
|
||||
calculate_deformation_gradient(ref->num_nbrs + 1, ref->points, res->mapping, scaled_points, ref->penrose, F, F_res);
|
||||
|
||||
if (P != NULL && U != NULL)
|
||||
polar_decomposition_3x3(F, false, U, P);
|
||||
}
|
||||
|
||||
if (mapping != NULL)
|
||||
for (int i=0;i<ref->num_nbrs + 1;i++)
|
||||
mapping[i] = ordering[res->mapping[i]];
|
||||
|
||||
double interatomic_distance = calculate_interatomic_distance(ref->type, res->scale);
|
||||
double lattice_constant = calculate_lattice_constant(ref->type, interatomic_distance);
|
||||
|
||||
if (p_interatomic_distance != NULL)
|
||||
*p_interatomic_distance = interatomic_distance;
|
||||
|
||||
if (p_lattice_constant != NULL)
|
||||
*p_lattice_constant = lattice_constant;
|
||||
|
||||
*p_rmsd = res->rmsd;
|
||||
*p_scale = res->scale;
|
||||
memcpy(q, res->q, 4 * sizeof(double));
|
||||
}
|
||||
|
||||
|
||||
extern bool ptm_initialized;
|
||||
|
||||
int ptm_index( ptm_local_handle_t local_handle, int32_t flags,
|
||||
int num_points, double (*unpermuted_points)[3], int32_t* unpermuted_numbers, bool topological_ordering,
|
||||
int32_t* p_type, int32_t* p_alloy_type, double* p_scale, double* p_rmsd, double* q, double* F, double* F_res,
|
||||
double* U, double* P, int8_t* mapping, double* p_interatomic_distance, double* p_lattice_constant)
|
||||
{
|
||||
assert(ptm_initialized);
|
||||
assert(num_points <= PTM_MAX_INPUT_POINTS);
|
||||
|
||||
if (flags & PTM_CHECK_SC)
|
||||
assert(num_points >= PTM_NUM_POINTS_SC);
|
||||
|
||||
if (flags & PTM_CHECK_BCC)
|
||||
assert(num_points >= PTM_NUM_POINTS_BCC);
|
||||
|
||||
if (flags & (PTM_CHECK_FCC | PTM_CHECK_HCP | PTM_CHECK_ICO))
|
||||
assert(num_points >= PTM_NUM_POINTS_FCC);
|
||||
|
||||
if (flags & (PTM_CHECK_DCUB | PTM_CHECK_DHEX))
|
||||
assert(num_points >= PTM_NUM_POINTS_DCUB);
|
||||
|
||||
int ret = 0;
|
||||
result_t res;
|
||||
res.ref_struct = NULL;
|
||||
res.rmsd = INFINITY;
|
||||
|
||||
int8_t ordering[PTM_MAX_INPUT_POINTS];
|
||||
double points[PTM_MAX_POINTS][3];
|
||||
int32_t numbers[PTM_MAX_POINTS];
|
||||
|
||||
int8_t dordering[PTM_MAX_INPUT_POINTS];
|
||||
double dpoints[PTM_MAX_POINTS][3];
|
||||
int32_t dnumbers[PTM_MAX_POINTS];
|
||||
|
||||
convexhull_t ch;
|
||||
double ch_points[PTM_MAX_INPUT_POINTS][3];
|
||||
|
||||
if (flags & (PTM_CHECK_SC | PTM_CHECK_FCC | PTM_CHECK_HCP | PTM_CHECK_ICO | PTM_CHECK_BCC))
|
||||
{
|
||||
int num_lpoints = std::min(std::min(PTM_MAX_POINTS, 20), num_points);
|
||||
order_points(local_handle, num_lpoints, unpermuted_points, unpermuted_numbers, topological_ordering, ordering, points, numbers);
|
||||
normalize_vertices(num_lpoints, points, ch_points);
|
||||
ch.ok = false;
|
||||
|
||||
if (flags & PTM_CHECK_SC)
|
||||
ret = match_general(&structure_sc, ch_points, points, &ch, &res);
|
||||
|
||||
if (flags & (PTM_CHECK_FCC | PTM_CHECK_HCP | PTM_CHECK_ICO))
|
||||
ret = match_fcc_hcp_ico(ch_points, points, flags, &ch, &res);
|
||||
|
||||
if (flags & PTM_CHECK_BCC)
|
||||
ret = match_general(&structure_bcc, ch_points, points, &ch, &res);
|
||||
}
|
||||
|
||||
if (flags & (PTM_CHECK_DCUB | PTM_CHECK_DHEX))
|
||||
{
|
||||
ret = calculate_diamond_neighbour_ordering(num_points, unpermuted_points, unpermuted_numbers, dordering, dpoints, dnumbers);
|
||||
if (ret == 0)
|
||||
{
|
||||
normalize_vertices(PTM_NUM_NBRS_DCUB + 1, dpoints, ch_points);
|
||||
ch.ok = false;
|
||||
|
||||
ret = match_dcub_dhex(ch_points, dpoints, flags, &ch, &res);
|
||||
}
|
||||
}
|
||||
|
||||
if (res.ref_struct != NULL && (res.ref_struct->type == PTM_MATCH_DCUB || res.ref_struct->type == PTM_MATCH_DHEX))
|
||||
{
|
||||
output_data( &res, num_points, unpermuted_numbers, dpoints, dnumbers, dordering,
|
||||
p_type, p_alloy_type, p_scale, p_rmsd, q, F, F_res,
|
||||
U, P, mapping, p_interatomic_distance, p_lattice_constant);
|
||||
}
|
||||
else
|
||||
{
|
||||
output_data( &res, num_points, unpermuted_numbers, points, numbers, ordering,
|
||||
p_type, p_alloy_type, p_scale, p_rmsd, q, F, F_res,
|
||||
U, P, mapping, p_interatomic_distance, p_lattice_constant);
|
||||
}
|
||||
|
||||
return PTM_NO_ERROR;
|
||||
}
|
||||
|
||||
71
src/PTM/ptm_initialize_data.cpp
Normal file
71
src/PTM/ptm_initialize_data.cpp
Normal file
@ -0,0 +1,71 @@
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <string.h>
|
||||
#include <cmath>
|
||||
#include <cfloat>
|
||||
#include <cassert>
|
||||
#include <algorithm>
|
||||
#include "ptm_initialize_data.h"
|
||||
|
||||
|
||||
static void make_facets_clockwise(int num_facets, int8_t (*facets)[3], const double (*points)[3])
|
||||
{
|
||||
double plane_normal[3];
|
||||
double origin[3] = {0, 0, 0};
|
||||
|
||||
for (int i = 0;i<num_facets;i++)
|
||||
add_facet(points, facets[i][0], facets[i][1], facets[i][2], facets[i], plane_normal, origin);
|
||||
}
|
||||
|
||||
static int initialize_graphs(const refdata_t* s, int8_t* colours)
|
||||
{
|
||||
for (int i = 0;i<s->num_graphs;i++)
|
||||
{
|
||||
int8_t code[2 * PTM_MAX_EDGES];
|
||||
int8_t degree[PTM_MAX_NBRS];
|
||||
int _max_degree = graph_degree(s->num_facets, s->graphs[i].facets, s->num_nbrs, degree);
|
||||
assert(_max_degree <= s->max_degree);
|
||||
|
||||
make_facets_clockwise(s->num_facets, s->graphs[i].facets, &s->points[1]);
|
||||
int ret = canonical_form_coloured(s->num_facets, s->graphs[i].facets, s->num_nbrs, degree, colours, s->graphs[i].canonical_labelling, (int8_t*)&code[0], &s->graphs[i].hash);
|
||||
if (ret != 0)
|
||||
return ret;
|
||||
}
|
||||
|
||||
return PTM_NO_ERROR;
|
||||
}
|
||||
|
||||
bool ptm_initialized = false;
|
||||
int ptm_initialize_global()
|
||||
{
|
||||
if (ptm_initialized)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
int8_t colours[PTM_MAX_POINTS] = {0};
|
||||
int8_t dcolours[PTM_MAX_POINTS] = {1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
|
||||
|
||||
int ret = initialize_graphs(&structure_sc, colours);
|
||||
ret |= initialize_graphs(&structure_fcc, colours);
|
||||
ret |= initialize_graphs(&structure_hcp, colours);
|
||||
ret |= initialize_graphs(&structure_ico, colours);
|
||||
ret |= initialize_graphs(&structure_bcc, colours);
|
||||
ret |= initialize_graphs(&structure_dcub, dcolours);
|
||||
ret |= initialize_graphs(&structure_dhex, dcolours);
|
||||
|
||||
if (ret == PTM_NO_ERROR)
|
||||
ptm_initialized = true;
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
ptm_local_handle_t ptm_initialize_local()
|
||||
{
|
||||
assert(ptm_initialized);
|
||||
return (ptm_local_handle_t)voronoi_initialize_local();
|
||||
}
|
||||
|
||||
void ptm_uninitialize_local(ptm_local_handle_t ptr)
|
||||
{
|
||||
voronoi_uninitialize_local(ptr);
|
||||
}
|
||||
|
||||
61
src/PTM/ptm_initialize_data.h
Normal file
61
src/PTM/ptm_initialize_data.h
Normal file
@ -0,0 +1,61 @@
|
||||
#ifndef PTM_INITIALIZE_DATA_H
|
||||
#define PTM_INITIALIZE_DATA_H
|
||||
|
||||
|
||||
#include "ptm_graph_data.h"
|
||||
#include "ptm_graph_tools.h"
|
||||
#include "ptm_deformation_gradient.h"
|
||||
#include "ptm_fundamental_mappings.h"
|
||||
#include "ptm_neighbour_ordering.h"
|
||||
#include "ptm_canonical_coloured.h"
|
||||
#include "ptm_convex_hull_incremental.h"
|
||||
|
||||
|
||||
typedef struct
|
||||
{
|
||||
int type;
|
||||
int num_nbrs;
|
||||
int num_facets;
|
||||
int max_degree;
|
||||
int num_graphs;
|
||||
int num_mappings;
|
||||
graph_t* graphs;
|
||||
const double (*points)[3];
|
||||
const double (*penrose)[3];
|
||||
const int8_t (*mapping)[PTM_MAX_POINTS];
|
||||
} refdata_t;
|
||||
|
||||
|
||||
//refdata_t structure_sc = { .type = PTM_MATCH_SC, .num_nbrs = 6, .num_facets = 8, .max_degree = 4, .num_graphs = NUM_SC_GRAPHS, .graphs = graphs_sc, .points = ptm_template_sc, .penrose = penrose_sc , .mapping = mapping_sc };
|
||||
const refdata_t structure_sc = { PTM_MATCH_SC, 6, 8, 4, NUM_SC_GRAPHS, NUM_CUBIC_MAPPINGS, graphs_sc, ptm_template_sc, penrose_sc, mapping_sc };
|
||||
const refdata_t structure_fcc = { PTM_MATCH_FCC, 12, 20, 6, NUM_FCC_GRAPHS, NUM_CUBIC_MAPPINGS, graphs_fcc, ptm_template_fcc, penrose_fcc, mapping_fcc };
|
||||
const refdata_t structure_hcp = { PTM_MATCH_HCP, 12, 20, 6, NUM_HCP_GRAPHS, NUM_HEX_MAPPINGS, graphs_hcp, ptm_template_hcp, penrose_hcp, mapping_hcp };
|
||||
const refdata_t structure_ico = { PTM_MATCH_ICO, 12, 20, 6, NUM_ICO_GRAPHS, NUM_ICO_MAPPINGS, graphs_ico, ptm_template_ico, penrose_ico, mapping_ico };
|
||||
const refdata_t structure_bcc = { PTM_MATCH_BCC, 14, 24, 8, NUM_BCC_GRAPHS, NUM_CUBIC_MAPPINGS, graphs_bcc, ptm_template_bcc, penrose_bcc, mapping_bcc };
|
||||
const refdata_t structure_dcub = { PTM_MATCH_DCUB, 16, 28, 8, NUM_DCUB_GRAPHS, NUM_DCUB_MAPPINGS, graphs_dcub, ptm_template_dcub, penrose_dcub, mapping_dcub };
|
||||
const refdata_t structure_dhex = { PTM_MATCH_DHEX, 16, 28, 8, NUM_DHEX_GRAPHS, NUM_DHEX_MAPPINGS, graphs_dhex, ptm_template_dhex, penrose_dhex, mapping_dhex };
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
typedef struct ptm_local_handle* ptm_local_handle_t;
|
||||
ptm_local_handle_t ptm_initialize_local();
|
||||
void ptm_uninitialize_local(ptm_local_handle_t ptr);
|
||||
|
||||
int ptm_initialize_global();
|
||||
|
||||
//------------------------------------
|
||||
// global initialization switch
|
||||
//------------------------------------
|
||||
extern bool ptm_initialized;
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
203
src/PTM/ptm_neighbour_ordering.cpp
Normal file
203
src/PTM/ptm_neighbour_ordering.cpp
Normal file
@ -0,0 +1,203 @@
|
||||
#include <cstdlib>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <cassert>
|
||||
#include <algorithm>
|
||||
#include "ptm_constants.h"
|
||||
#include "ptm_voronoi_cell.h"
|
||||
using namespace voro;
|
||||
|
||||
|
||||
|
||||
typedef struct
|
||||
{
|
||||
double area;
|
||||
double dist;
|
||||
int index;
|
||||
} sorthelper_t;
|
||||
|
||||
static bool sorthelper_compare(sorthelper_t const& a, sorthelper_t const& b)
|
||||
{
|
||||
if (a.area > b.area)
|
||||
return true;
|
||||
|
||||
if (a.area < b.area)
|
||||
return false;
|
||||
|
||||
if (a.dist < b.dist)
|
||||
return true;
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
//todo: change voronoi code to return errors rather than exiting
|
||||
static int calculate_voronoi_face_areas(int num_points, const double (*_points)[3], double* normsq, double max_norm, voronoicell_neighbor* v, std::vector<int>& nbr_indices, std::vector<double>& face_areas)
|
||||
{
|
||||
const double k = 1000 * max_norm; //todo: reduce this constant
|
||||
v->init(-k,k,-k,k,-k,k);
|
||||
|
||||
for (int i=1;i<num_points;i++)
|
||||
{
|
||||
double x = _points[i][0] - _points[0][0];
|
||||
double y = _points[i][1] - _points[0][1];
|
||||
double z = _points[i][2] - _points[0][2];
|
||||
v->nplane(x,y,z,normsq[i],i);
|
||||
}
|
||||
|
||||
v->neighbors(nbr_indices);
|
||||
v->face_areas(face_areas);
|
||||
return 0;
|
||||
}
|
||||
|
||||
int calculate_neighbour_ordering(void* _voronoi_handle, int num_points, const double (*_points)[3], int8_t* ordering)
|
||||
{
|
||||
assert(num_points <= PTM_MAX_INPUT_POINTS);
|
||||
|
||||
voronoicell_neighbor* voronoi_handle = (voronoicell_neighbor*)_voronoi_handle;
|
||||
|
||||
double max_norm = 0;
|
||||
double points[PTM_MAX_INPUT_POINTS][3];
|
||||
double normsq[PTM_MAX_INPUT_POINTS];
|
||||
for (int i = 0;i<num_points;i++)
|
||||
{
|
||||
double x = _points[i][0] - _points[0][0];
|
||||
double y = _points[i][1] - _points[0][1];
|
||||
double z = _points[i][2] - _points[0][2];
|
||||
points[i][0] = x;
|
||||
points[i][1] = y;
|
||||
points[i][2] = z;
|
||||
|
||||
normsq[i] = x*x + y*y + z*z;
|
||||
max_norm = std::max(max_norm, normsq[i]);
|
||||
#ifdef DEBUG
|
||||
printf("point %d: %f\t%f\t%f\t%f\n", i, x, y, z, x*x + y*y + z*z);
|
||||
#endif
|
||||
}
|
||||
|
||||
max_norm = sqrt(max_norm);
|
||||
|
||||
std::vector<int> nbr_indices(num_points + 6);
|
||||
std::vector<double> face_areas(num_points + 6);
|
||||
int ret = calculate_voronoi_face_areas(num_points, points, normsq, max_norm, voronoi_handle, nbr_indices, face_areas);
|
||||
if (ret != 0)
|
||||
return ret;
|
||||
|
||||
double areas[PTM_MAX_INPUT_POINTS];
|
||||
memset(areas, 0, num_points * sizeof(double));
|
||||
areas[0] = INFINITY;
|
||||
for (size_t i=0;i<nbr_indices.size();i++)
|
||||
{
|
||||
int index = nbr_indices[i];
|
||||
if (index > 0)
|
||||
areas[index] = face_areas[i];
|
||||
}
|
||||
|
||||
sorthelper_t data[PTM_MAX_INPUT_POINTS];
|
||||
for (int i=0;i<num_points;i++)
|
||||
{
|
||||
assert(areas[i] == areas[i]);
|
||||
data[i].area = areas[i];
|
||||
data[i].dist = normsq[i];
|
||||
data[i].index = i;
|
||||
}
|
||||
|
||||
std::sort(data, data + num_points, &sorthelper_compare);
|
||||
|
||||
#ifdef DEBUG
|
||||
for (int i=0;i<num_points;i++)
|
||||
printf("%d %f\n", data[i].index, data[i].area);
|
||||
#endif
|
||||
|
||||
for (int i=0;i<num_points;i++)
|
||||
ordering[i] = data[i].index;
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
void* voronoi_initialize_local()
|
||||
{
|
||||
voronoicell_neighbor* ptr = new voronoicell_neighbor;
|
||||
return (void*)ptr;
|
||||
}
|
||||
|
||||
void voronoi_uninitialize_local(void* _ptr)
|
||||
{
|
||||
voronoicell_neighbor* ptr = (voronoicell_neighbor*)_ptr;
|
||||
delete ptr;
|
||||
}
|
||||
|
||||
|
||||
typedef struct
|
||||
{
|
||||
double dist;
|
||||
int p;
|
||||
int index;
|
||||
} diamond_t;
|
||||
|
||||
static bool diamond_compare(diamond_t const& a, diamond_t const& b)
|
||||
{
|
||||
return a.dist < b.dist;
|
||||
}
|
||||
|
||||
int calculate_diamond_neighbour_ordering( int num_points, double (*unpermuted_points)[3], int32_t* unpermuted_numbers,
|
||||
int8_t* ordering, double (*points)[3], int32_t* numbers)
|
||||
{
|
||||
assert(num_points <= PTM_MAX_INPUT_POINTS);
|
||||
|
||||
diamond_t data[4 * (PTM_MAX_INPUT_POINTS - 5)];
|
||||
int index = 0;
|
||||
for (int i=5;i<num_points;i++)
|
||||
{
|
||||
for (int j=1;j<5;j++)
|
||||
{
|
||||
double dx = unpermuted_points[i][0] - unpermuted_points[j][0];
|
||||
double dy = unpermuted_points[i][1] - unpermuted_points[j][1];
|
||||
double dz = unpermuted_points[i][2] - unpermuted_points[j][2];
|
||||
|
||||
double d = dx*dx + dy*dy + dz*dz;
|
||||
|
||||
data[index].p = j - 1;
|
||||
data[index].index = i;
|
||||
data[index].dist = d;
|
||||
index++;
|
||||
}
|
||||
}
|
||||
int n = index;
|
||||
|
||||
std::sort(data, data + n, &diamond_compare);
|
||||
|
||||
for (index=0;index<5;index++)
|
||||
ordering[index] = index;
|
||||
|
||||
int num_found = 0;
|
||||
bool hit[PTM_MAX_INPUT_POINTS] = {0};
|
||||
int counts[4] = {0};
|
||||
for (int i=0;i<n;i++)
|
||||
{
|
||||
int p = data[i].p;
|
||||
int q = data[i].index;
|
||||
if (hit[q] || counts[p] >= 3)
|
||||
continue;
|
||||
|
||||
ordering[1 + 4 + 3 * p + counts[p]] = q;
|
||||
counts[p]++;
|
||||
index++;
|
||||
num_found++;
|
||||
if (num_found >= 12)
|
||||
break;
|
||||
}
|
||||
|
||||
if (num_found != 12)
|
||||
return -1;
|
||||
|
||||
for (int i=0;i<PTM_NUM_NBRS_DCUB+1;i++)
|
||||
{
|
||||
memcpy(points[i], &unpermuted_points[ordering[i]], 3 * sizeof(double));
|
||||
|
||||
if (unpermuted_numbers != NULL)
|
||||
numbers[i] = unpermuted_numbers[ordering[i]];
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
13
src/PTM/ptm_neighbour_ordering.h
Normal file
13
src/PTM/ptm_neighbour_ordering.h
Normal file
@ -0,0 +1,13 @@
|
||||
#ifndef PTM_NEIGHBOUR_ORDERING_H
|
||||
#define PTM_NEIGHBOUR_ORDERING_H
|
||||
|
||||
int calculate_neighbour_ordering(void* voronoi_handle, int num_points, const double (*_points)[3], int8_t* ordering);
|
||||
|
||||
int calculate_diamond_neighbour_ordering( int num_points, double (*unpermuted_points)[3], int32_t* unpermuted_numbers,
|
||||
int8_t* ordering, double (*points)[3], int32_t* numbers);
|
||||
|
||||
void* voronoi_initialize_local();
|
||||
void voronoi_uninitialize_local(void* ptr);
|
||||
|
||||
#endif
|
||||
|
||||
55
src/PTM/ptm_normalize_vertices.cpp
Normal file
55
src/PTM/ptm_normalize_vertices.cpp
Normal file
@ -0,0 +1,55 @@
|
||||
#include <cmath>
|
||||
|
||||
|
||||
void subtract_barycentre(int num, double (*points)[3], double (*normalized)[3])
|
||||
{
|
||||
//calculate barycentre
|
||||
double sum[3] = {0, 0, 0};
|
||||
for (int i=0;i<num;i++)
|
||||
{
|
||||
sum[0] += points[i][0];
|
||||
sum[1] += points[i][1];
|
||||
sum[2] += points[i][2];
|
||||
}
|
||||
|
||||
sum[0] /= num;
|
||||
sum[1] /= num;
|
||||
sum[2] /= num;
|
||||
|
||||
//subtract barycentre
|
||||
for (int i=0;i<num;i++)
|
||||
{
|
||||
normalized[i][0] = points[i][0] - sum[0];
|
||||
normalized[i][1] = points[i][1] - sum[1];
|
||||
normalized[i][2] = points[i][2] - sum[2];
|
||||
}
|
||||
}
|
||||
|
||||
double normalize_vertices(int num, double (*points)[3], double (*normalized)[3])
|
||||
{
|
||||
subtract_barycentre(num, points, normalized);
|
||||
|
||||
//calculate mean length
|
||||
double scale = 0.0;
|
||||
for (int i=1;i<num;i++)
|
||||
{
|
||||
double x = normalized[i][0];
|
||||
double y = normalized[i][1];
|
||||
double z = normalized[i][2];
|
||||
|
||||
double norm = sqrt(x*x + y*y + z*z);
|
||||
scale += norm;
|
||||
}
|
||||
scale /= num;
|
||||
|
||||
//scale vertices such that mean length is 1
|
||||
for (int i=0;i<num;i++)
|
||||
{
|
||||
normalized[i][0] /= scale;
|
||||
normalized[i][1] /= scale;
|
||||
normalized[i][2] /= scale;
|
||||
}
|
||||
|
||||
return scale;
|
||||
}
|
||||
|
||||
8
src/PTM/ptm_normalize_vertices.h
Normal file
8
src/PTM/ptm_normalize_vertices.h
Normal file
@ -0,0 +1,8 @@
|
||||
#ifndef PTM_NORMALIZE_VERTICES_H
|
||||
#define PTM_NORMALIZE_VERTICES_H
|
||||
|
||||
void subtract_barycentre(int num, double (*points)[3], double (*normalized)[3]);
|
||||
double normalize_vertices(int num, double (*points)[3], double (*normalized)[3]);
|
||||
|
||||
#endif
|
||||
|
||||
337
src/PTM/ptm_polar.cpp
Normal file
337
src/PTM/ptm_polar.cpp
Normal file
@ -0,0 +1,337 @@
|
||||
/*******************************************************************************
|
||||
* -/_|:|_|_\-
|
||||
*
|
||||
* This code is a modification of D.L. Theobald's QCP rotation code.
|
||||
* It has been adapted to calculate the polar decomposition of a 3x3 matrix
|
||||
* Adaption by P.M. Larsen
|
||||
*
|
||||
* Original Author(s): Douglas L. Theobald
|
||||
* Department of Biochemistry
|
||||
* MS 009
|
||||
* Brandeis University
|
||||
* 415 South St
|
||||
* Waltham, MA 02453
|
||||
* USA
|
||||
*
|
||||
* dtheobald@brandeis.edu
|
||||
*
|
||||
* Pu Liu
|
||||
* Johnson & Johnson Pharmaceutical Research and Development, L.L.C.
|
||||
* 665 Stockton Drive
|
||||
* Exton, PA 19341
|
||||
* USA
|
||||
*
|
||||
* pliu24@its.jnj.com
|
||||
*
|
||||
*
|
||||
* If you use this QCP rotation calculation method in a publication, please
|
||||
* reference:
|
||||
*
|
||||
* Douglas L. Theobald (2005)
|
||||
* "Rapid calculation of RMSD using a quaternion-based characteristic
|
||||
* polynomial."
|
||||
* Acta Crystallographica A 61(4):478-480.
|
||||
*
|
||||
* Pu Liu, Dmitris K. Agrafiotis, and Douglas L. Theobald (2009)
|
||||
* "Fast determination of the optimal rotational matrix for macromolecular
|
||||
* superpositions."
|
||||
* Journal of Computational Chemistry 31(7):1561-1563.
|
||||
*
|
||||
*
|
||||
* Copyright (c) 2009-2013 Pu Liu and Douglas L. Theobald
|
||||
* All rights reserved.
|
||||
*
|
||||
* Redistribution and use in source and binary forms, with or without modification, are permitted
|
||||
* provided that the following conditions are met:
|
||||
*
|
||||
* * Redistributions of source code must retain the above copyright notice, this list of
|
||||
* conditions and the following disclaimer.
|
||||
* * Redistributions in binary form must reproduce the above copyright notice, this list
|
||||
* of conditions and the following disclaimer in the documentation and/or other materials
|
||||
* provided with the distribution.
|
||||
* * Neither the name of the <ORGANIZATION> nor the names of its contributors may be used to
|
||||
* endorse or promote products derived from this software without specific prior written
|
||||
* permission.
|
||||
*
|
||||
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
||||
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
||||
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
|
||||
* PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
|
||||
* HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
|
||||
* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
|
||||
* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
||||
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
||||
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
||||
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
|
||||
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
||||
*
|
||||
* Source: started anew.
|
||||
*
|
||||
* Change History:
|
||||
* 2009/04/13 Started source
|
||||
* 2010/03/28 Modified FastCalcRMSDAndRotation() to handle tiny qsqr
|
||||
* If trying all rows of the adjoint still gives too small
|
||||
* qsqr, then just return identity matrix. (DLT)
|
||||
* 2010/06/30 Fixed prob in assigning A[9] = 0 in InnerProduct()
|
||||
* invalid mem access
|
||||
* 2011/02/21 Made CenterCoords use weights
|
||||
* 2011/05/02 Finally changed CenterCoords declaration in qcprot.h
|
||||
* Also changed some functions to static
|
||||
* 2011/07/08 put in fabs() to fix taking sqrt of small neg numbers, fp error
|
||||
* 2012/07/26 minor changes to comments and main.c, more info (v.1.4)
|
||||
*
|
||||
* 2016/05/29 QCP method adapted for polar decomposition of a 3x3 matrix,
|
||||
* for use in Polyhedral Template Matching.
|
||||
*
|
||||
******************************************************************************/
|
||||
|
||||
#include <cmath>
|
||||
#include <algorithm>
|
||||
#include <string.h>
|
||||
#include "ptm_quat.h"
|
||||
|
||||
|
||||
static void matmul_3x3(double* A, double* x, double* b)
|
||||
{
|
||||
b[0] = A[0] * x[0] + A[1] * x[3] + A[2] * x[6];
|
||||
b[3] = A[3] * x[0] + A[4] * x[3] + A[5] * x[6];
|
||||
b[6] = A[6] * x[0] + A[7] * x[3] + A[8] * x[6];
|
||||
|
||||
b[1] = A[0] * x[1] + A[1] * x[4] + A[2] * x[7];
|
||||
b[4] = A[3] * x[1] + A[4] * x[4] + A[5] * x[7];
|
||||
b[7] = A[6] * x[1] + A[7] * x[4] + A[8] * x[7];
|
||||
|
||||
b[2] = A[0] * x[2] + A[1] * x[5] + A[2] * x[8];
|
||||
b[5] = A[3] * x[2] + A[4] * x[5] + A[5] * x[8];
|
||||
b[8] = A[6] * x[2] + A[7] * x[5] + A[8] * x[8];
|
||||
}
|
||||
|
||||
static double matrix_determinant_3x3(double* A)
|
||||
{
|
||||
return A[0] * (A[4]*A[8] - A[5]*A[7])
|
||||
- A[1] * (A[3]*A[8] - A[5]*A[6])
|
||||
+ A[2] * (A[3]*A[7] - A[4]*A[6]);
|
||||
}
|
||||
|
||||
static void flip_matrix(double* A)
|
||||
{
|
||||
for (int i=0;i<9;i++)
|
||||
A[i] = -A[i];
|
||||
}
|
||||
|
||||
static bool optimal_quaternion(double* A, bool polar, double E0, double* p_nrmsdsq, double* qopt)
|
||||
{
|
||||
const double evecprec = 1e-6;
|
||||
const double evalprec = 1e-11;
|
||||
|
||||
double Sxx = A[0], Sxy = A[1], Sxz = A[2],
|
||||
Syx = A[3], Syy = A[4], Syz = A[5],
|
||||
Szx = A[6], Szy = A[7], Szz = A[8];
|
||||
|
||||
double Sxx2 = Sxx * Sxx, Syy2 = Syy * Syy, Szz2 = Szz * Szz,
|
||||
Sxy2 = Sxy * Sxy, Syz2 = Syz * Syz, Sxz2 = Sxz * Sxz,
|
||||
Syx2 = Syx * Syx, Szy2 = Szy * Szy, Szx2 = Szx * Szx;
|
||||
|
||||
double fnorm_squared = Sxx2 + Syy2 + Szz2 + Sxy2 + Syz2 + Sxz2 + Syx2 + Szy2 + Szx2;
|
||||
|
||||
double SyzSzymSyySzz2 = 2.0 * (Syz * Szy - Syy * Szz);
|
||||
double Sxx2Syy2Szz2Syz2Szy2 = Syy2 + Szz2 - Sxx2 + Syz2 + Szy2;
|
||||
double SxzpSzx = Sxz + Szx;
|
||||
double SyzpSzy = Syz + Szy;
|
||||
double SxypSyx = Sxy + Syx;
|
||||
double SyzmSzy = Syz - Szy;
|
||||
double SxzmSzx = Sxz - Szx;
|
||||
double SxymSyx = Sxy - Syx;
|
||||
double SxxpSyy = Sxx + Syy;
|
||||
double SxxmSyy = Sxx - Syy;
|
||||
double Sxy2Sxz2Syx2Szx2 = Sxy2 + Sxz2 - Syx2 - Szx2;
|
||||
|
||||
double C[3];
|
||||
C[0] = Sxy2Sxz2Syx2Szx2 * Sxy2Sxz2Syx2Szx2
|
||||
+ (Sxx2Syy2Szz2Syz2Szy2 + SyzSzymSyySzz2) * (Sxx2Syy2Szz2Syz2Szy2 - SyzSzymSyySzz2)
|
||||
+ (-(SxzpSzx)*(SyzmSzy)+(SxymSyx)*(SxxmSyy-Szz)) * (-(SxzmSzx)*(SyzpSzy)+(SxymSyx)*(SxxmSyy+Szz))
|
||||
+ (-(SxzpSzx)*(SyzpSzy)-(SxypSyx)*(SxxpSyy-Szz)) * (-(SxzmSzx)*(SyzmSzy)-(SxypSyx)*(SxxpSyy+Szz))
|
||||
+ (+(SxypSyx)*(SyzpSzy)+(SxzpSzx)*(SxxmSyy+Szz)) * (-(SxymSyx)*(SyzmSzy)+(SxzpSzx)*(SxxpSyy+Szz))
|
||||
+ (+(SxypSyx)*(SyzmSzy)+(SxzmSzx)*(SxxmSyy-Szz)) * (-(SxymSyx)*(SyzpSzy)+(SxzmSzx)*(SxxpSyy-Szz));
|
||||
|
||||
C[1] = 8.0 * (Sxx*Syz*Szy + Syy*Szx*Sxz + Szz*Sxy*Syx - Sxx*Syy*Szz - Syz*Szx*Sxy - Szy*Syx*Sxz);
|
||||
C[2] = -2.0 * fnorm_squared;
|
||||
|
||||
//Newton-Raphson
|
||||
double mxEigenV = polar ? sqrt(3 * fnorm_squared) : E0;
|
||||
if (mxEigenV > evalprec)
|
||||
{
|
||||
for (int i=0;i<50;i++)
|
||||
{
|
||||
double oldg = mxEigenV;
|
||||
double x2 = mxEigenV*mxEigenV;
|
||||
double b = (x2 + C[2])*mxEigenV;
|
||||
double a = b + C[1];
|
||||
double delta = ((a * mxEigenV + C[0]) / (2 * x2 * mxEigenV + b + a));
|
||||
mxEigenV -= delta;
|
||||
if (fabs(mxEigenV - oldg) < fabs(evalprec * mxEigenV))
|
||||
break;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
mxEigenV = 0.0;
|
||||
}
|
||||
|
||||
(*p_nrmsdsq) = std::max(0.0, 2.0 * (E0 - mxEigenV));
|
||||
|
||||
double a11 = SxxpSyy + Szz - mxEigenV;
|
||||
double a12 = SyzmSzy;
|
||||
double a13 = -SxzmSzx;
|
||||
double a14 = SxymSyx;
|
||||
|
||||
double a21 = SyzmSzy;
|
||||
double a22 = SxxmSyy - Szz -mxEigenV;
|
||||
double a23 = SxypSyx;
|
||||
double a24 = SxzpSzx;
|
||||
|
||||
double a31 = a13;
|
||||
double a32 = a23;
|
||||
double a33 = Syy - Sxx - Szz - mxEigenV;
|
||||
double a34 = SyzpSzy;
|
||||
|
||||
double a41 = a14;
|
||||
double a42 = a24;
|
||||
double a43 = a34;
|
||||
double a44 = Szz - SxxpSyy - mxEigenV;
|
||||
|
||||
double a3344_4334 = a33 * a44 - a43 * a34;
|
||||
double a3244_4234 = a32 * a44 - a42 * a34;
|
||||
double a3243_4233 = a32 * a43 - a42 * a33;
|
||||
double a3143_4133 = a31 * a43 - a41 * a33;
|
||||
double a3144_4134 = a31 * a44 - a41 * a34;
|
||||
double a3142_4132 = a31 * a42 - a41 * a32;
|
||||
double a1324_1423 = a13 * a24 - a14 * a23;
|
||||
double a1224_1422 = a12 * a24 - a14 * a22;
|
||||
double a1223_1322 = a12 * a23 - a13 * a22;
|
||||
double a1124_1421 = a11 * a24 - a14 * a21;
|
||||
double a1123_1321 = a11 * a23 - a13 * a21;
|
||||
double a1122_1221 = a11 * a22 - a12 * a21;
|
||||
|
||||
double q[4][4];
|
||||
q[0][0] = a12 * a3344_4334 - a13 * a3244_4234 + a14 * a3243_4233;
|
||||
q[0][1] = -a11 * a3344_4334 + a13 * a3144_4134 - a14 * a3143_4133;
|
||||
q[0][2] = a11 * a3244_4234 - a12 * a3144_4134 + a14 * a3142_4132;
|
||||
q[0][3] = -a11 * a3243_4233 + a12 * a3143_4133 - a13 * a3142_4132;
|
||||
|
||||
q[1][0] = a22 * a3344_4334 - a23 * a3244_4234 + a24 * a3243_4233;
|
||||
q[1][1] = -a21 * a3344_4334 + a23 * a3144_4134 - a24 * a3143_4133;
|
||||
q[1][2] = a21 * a3244_4234 - a22 * a3144_4134 + a24 * a3142_4132;
|
||||
q[1][3] = -a21 * a3243_4233 + a22 * a3143_4133 - a23 * a3142_4132;
|
||||
|
||||
q[2][0] = a32 * a1324_1423 - a33 * a1224_1422 + a34 * a1223_1322;
|
||||
q[2][1] = -a31 * a1324_1423 + a33 * a1124_1421 - a34 * a1123_1321;
|
||||
q[2][2] = a31 * a1224_1422 - a32 * a1124_1421 + a34 * a1122_1221;
|
||||
q[2][3] = -a31 * a1223_1322 + a32 * a1123_1321 - a33 * a1122_1221;
|
||||
|
||||
q[3][0] = a42 * a1324_1423 - a43 * a1224_1422 + a44 * a1223_1322;
|
||||
q[3][1] = -a41 * a1324_1423 + a43 * a1124_1421 - a44 * a1123_1321;
|
||||
q[3][2] = a41 * a1224_1422 - a42 * a1124_1421 + a44 * a1122_1221;
|
||||
q[3][3] = -a41 * a1223_1322 + a42 * a1123_1321 - a43 * a1122_1221;
|
||||
|
||||
double qsqr[4];
|
||||
for (int i=0;i<4;i++)
|
||||
qsqr[i] = q[i][0]*q[i][0] + q[i][1]*q[i][1] + q[i][2]*q[i][2] + q[i][3]*q[i][3];
|
||||
|
||||
int bi = 0;
|
||||
double max = 0;
|
||||
for (int i=0;i<4;i++)
|
||||
{
|
||||
if (qsqr[i] > max)
|
||||
{
|
||||
bi = i;
|
||||
max = qsqr[i];
|
||||
}
|
||||
}
|
||||
|
||||
bool too_small = false;
|
||||
if (qsqr[bi] < evecprec)
|
||||
{
|
||||
//if qsqr is still too small, return the identity rotation.
|
||||
q[bi][0] = 1;
|
||||
q[bi][1] = 0;
|
||||
q[bi][2] = 0;
|
||||
q[bi][3] = 0;
|
||||
too_small = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
double normq = sqrt(qsqr[bi]);
|
||||
q[bi][0] /= normq;
|
||||
q[bi][1] /= normq;
|
||||
q[bi][2] /= normq;
|
||||
q[bi][3] /= normq;
|
||||
}
|
||||
|
||||
memcpy(qopt, q[bi], 4 * sizeof(double));
|
||||
return !too_small;
|
||||
}
|
||||
|
||||
int polar_decomposition_3x3(double* _A, bool right_sided, double* U, double* P)
|
||||
{
|
||||
double A[9];
|
||||
memcpy(A, _A, 9 * sizeof(double));
|
||||
|
||||
double det = matrix_determinant_3x3(A);
|
||||
if (det < 0)
|
||||
flip_matrix(A);
|
||||
|
||||
double q[4];
|
||||
double nrmsdsq = 0;
|
||||
optimal_quaternion(A, true, -1, &nrmsdsq, q);
|
||||
q[0] = -q[0];
|
||||
quaternion_to_rotation_matrix(q, U);
|
||||
|
||||
if (det < 0)
|
||||
flip_matrix(U);
|
||||
|
||||
double UT[9] = {U[0], U[3], U[6], U[1], U[4], U[7], U[2], U[5], U[8]};
|
||||
|
||||
if (right_sided)
|
||||
matmul_3x3(UT, _A, P);
|
||||
else
|
||||
matmul_3x3(_A, UT, P);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
void InnerProduct(double *A, int num, const double (*coords1)[3], double (*coords2)[3], int8_t* permutation)
|
||||
{
|
||||
A[0] = A[1] = A[2] = A[3] = A[4] = A[5] = A[6] = A[7] = A[8] = 0.0;
|
||||
|
||||
for (int i = 0; i < num; ++i)
|
||||
{
|
||||
double x1 = coords1[i][0];
|
||||
double y1 = coords1[i][1];
|
||||
double z1 = coords1[i][2];
|
||||
|
||||
double x2 = coords2[permutation[i]][0];
|
||||
double y2 = coords2[permutation[i]][1];
|
||||
double z2 = coords2[permutation[i]][2];
|
||||
|
||||
A[0] += x1 * x2;
|
||||
A[1] += x1 * y2;
|
||||
A[2] += x1 * z2;
|
||||
|
||||
A[3] += y1 * x2;
|
||||
A[4] += y1 * y2;
|
||||
A[5] += y1 * z2;
|
||||
|
||||
A[6] += z1 * x2;
|
||||
A[7] += z1 * y2;
|
||||
A[8] += z1 * z2;
|
||||
}
|
||||
}
|
||||
|
||||
int FastCalcRMSDAndRotation(double *A, double E0, double *p_nrmsdsq, double *q, double* U)
|
||||
{
|
||||
optimal_quaternion(A, false, E0, p_nrmsdsq, q);
|
||||
quaternion_to_rotation_matrix(q, U);
|
||||
return 0;
|
||||
}
|
||||
|
||||
12
src/PTM/ptm_polar.h
Normal file
12
src/PTM/ptm_polar.h
Normal file
@ -0,0 +1,12 @@
|
||||
#ifndef PTM_POLAR_H
|
||||
#define PTM_POLAR_H
|
||||
|
||||
#include <stdint.h>
|
||||
#include <stdbool.h>
|
||||
|
||||
int polar_decomposition_3x3(double* _A, bool right_sided, double* U, double* P);
|
||||
void InnerProduct(double *A, int num, const double (*coords1)[3], double (*coords2)[3], int8_t* permutation);
|
||||
int FastCalcRMSDAndRotation(double *A, double E0, double *p_nrmsdsq, double *q, double* U);
|
||||
|
||||
#endif
|
||||
|
||||
396
src/PTM/ptm_quat.cpp
Normal file
396
src/PTM/ptm_quat.cpp
Normal file
@ -0,0 +1,396 @@
|
||||
#include <string.h>
|
||||
#include <cmath>
|
||||
#include <cfloat>
|
||||
|
||||
|
||||
#define SIGN(x) (x >= 0 ? 1 : -1)
|
||||
#define MIN(X, Y) (((X) < (Y)) ? (X) : (Y))
|
||||
#define MAX(X, Y) (((X) > (Y)) ? (X) : (Y))
|
||||
|
||||
|
||||
#define SQRT_2 1.4142135623730951454746218587388284504414
|
||||
#define HALF_SQRT_2 0.7071067811865474617150084668537601828575
|
||||
|
||||
#define PHI 1.6180339887498949025257388711906969547272
|
||||
#define HALF_PHI 0.8090169943749474512628694355953484773636
|
||||
|
||||
#define INV_PHI 0.6180339887498947915034364086750429123640
|
||||
#define HALF_INV_PHI 0.3090169943749473957517182043375214561820
|
||||
|
||||
#define SQRT_5_ 2.23606797749978969640917366873127623544061835961152572427089
|
||||
#define SQRT_2_3 0.8164965809277260344600790631375275552273
|
||||
#define SQRT_1_6 0.4082482904638630172300395315687637776136
|
||||
|
||||
|
||||
double generator_cubic[24][4] = { {1, 0, 0, 0 },
|
||||
{0, 1, 0, 0 },
|
||||
{0, 0, 1, 0 },
|
||||
{0, 0, 0, 1 },
|
||||
{0.5, 0.5, 0.5, 0.5 },
|
||||
{0.5, 0.5, -0.5, 0.5 },
|
||||
{0.5, -0.5, 0.5, 0.5 },
|
||||
{0.5, -0.5, -0.5, 0.5 },
|
||||
{-0.5, 0.5, 0.5, 0.5 },
|
||||
{-0.5, 0.5, -0.5, 0.5 },
|
||||
{-0.5, -0.5, 0.5, 0.5 },
|
||||
{-0.5, -0.5, -0.5, 0.5 },
|
||||
{HALF_SQRT_2, HALF_SQRT_2, 0, 0 },
|
||||
{HALF_SQRT_2, 0, HALF_SQRT_2, 0 },
|
||||
{HALF_SQRT_2, 0, 0, HALF_SQRT_2 },
|
||||
{-HALF_SQRT_2, HALF_SQRT_2, 0, 0 },
|
||||
{-HALF_SQRT_2, 0, HALF_SQRT_2, 0 },
|
||||
{-HALF_SQRT_2, 0, 0, HALF_SQRT_2 },
|
||||
{0, HALF_SQRT_2, HALF_SQRT_2, 0 },
|
||||
{0, HALF_SQRT_2, 0, HALF_SQRT_2 },
|
||||
{0, 0, HALF_SQRT_2, HALF_SQRT_2 },
|
||||
{0, -HALF_SQRT_2, HALF_SQRT_2, 0 },
|
||||
{0, -HALF_SQRT_2, 0, HALF_SQRT_2 },
|
||||
{0, 0, -HALF_SQRT_2, HALF_SQRT_2 } };
|
||||
|
||||
double generator_diamond_cubic[12][4] = { {1, 0, 0, 0 },
|
||||
{0, 1, 0, 0 },
|
||||
{0, 0, 1, 0 },
|
||||
{0, 0, 0, 1 },
|
||||
{0.5, 0.5, 0.5, 0.5 },
|
||||
{0.5, 0.5, -0.5, 0.5 },
|
||||
{0.5, -0.5, 0.5, 0.5 },
|
||||
{0.5, -0.5, -0.5, 0.5 },
|
||||
{-0.5, 0.5, 0.5, 0.5 },
|
||||
{-0.5, 0.5, -0.5, 0.5 },
|
||||
{-0.5, -0.5, 0.5, 0.5 },
|
||||
{-0.5, -0.5, -0.5, 0.5 } };
|
||||
|
||||
double generator_hcp[6][4] = { {1, 0, 0, 0},
|
||||
{0.5, 0.5, 0.5, 0.5},
|
||||
{0.5, -0.5, -0.5, -0.5},
|
||||
{0, SQRT_2_3, -SQRT_1_6, -SQRT_1_6},
|
||||
{0, SQRT_1_6, -SQRT_2_3, SQRT_1_6},
|
||||
{0, SQRT_1_6, SQRT_1_6, -SQRT_2_3} };
|
||||
|
||||
double generator_diamond_hexagonal[3][4] = { {1, 0, 0, 0},
|
||||
{0.5, 0.5, 0.5, 0.5},
|
||||
{0.5, -0.5, -0.5, -0.5} };
|
||||
|
||||
double generator_icosahedral[60][4] = { {1, 0, 0, 0},
|
||||
{HALF_PHI, -HALF_INV_PHI, -0.5, 0},
|
||||
{HALF_PHI, 0, -HALF_INV_PHI, -0.5},
|
||||
{HALF_PHI, -0.5, 0, -HALF_INV_PHI},
|
||||
{HALF_PHI, HALF_INV_PHI, -0.5, 0},
|
||||
{HALF_PHI, 0, HALF_INV_PHI, -0.5},
|
||||
{HALF_PHI, -0.5, 0, HALF_INV_PHI},
|
||||
{HALF_PHI, 0.5, 0, -HALF_INV_PHI},
|
||||
{HALF_PHI, 0, -HALF_INV_PHI, 0.5},
|
||||
{HALF_PHI, -HALF_INV_PHI, 0.5, 0},
|
||||
{HALF_PHI, 0, HALF_INV_PHI, 0.5},
|
||||
{HALF_PHI, HALF_INV_PHI, 0.5, 0},
|
||||
{HALF_PHI, 0.5, 0, HALF_INV_PHI},
|
||||
{0.5, HALF_PHI, -HALF_INV_PHI, 0},
|
||||
{0.5, HALF_PHI, HALF_INV_PHI, 0},
|
||||
{0.5, 0.5, 0.5, 0.5},
|
||||
{0.5, 0.5, 0.5, -0.5},
|
||||
{0.5, 0.5, -0.5, 0.5},
|
||||
{0.5, 0.5, -0.5, -0.5},
|
||||
{0.5, HALF_INV_PHI, 0, HALF_PHI},
|
||||
{0.5, HALF_INV_PHI, 0, -HALF_PHI},
|
||||
{0.5, 0, HALF_PHI, -HALF_INV_PHI},
|
||||
{0.5, 0, HALF_PHI, HALF_INV_PHI},
|
||||
{0.5, 0, -HALF_PHI, -HALF_INV_PHI},
|
||||
{0.5, 0, -HALF_PHI, HALF_INV_PHI},
|
||||
{0.5, -HALF_INV_PHI, 0, HALF_PHI},
|
||||
{0.5, -HALF_INV_PHI, 0, -HALF_PHI},
|
||||
{0.5, -0.5, 0.5, 0.5},
|
||||
{0.5, -0.5, 0.5, -0.5},
|
||||
{0.5, -0.5, -0.5, 0.5},
|
||||
{0.5, -0.5, -0.5, -0.5},
|
||||
{0.5, -HALF_PHI, -HALF_INV_PHI, 0},
|
||||
{0.5, -HALF_PHI, HALF_INV_PHI, 0},
|
||||
{HALF_INV_PHI, -HALF_PHI, 0, -0.5},
|
||||
{HALF_INV_PHI, 0, -0.5, -HALF_PHI},
|
||||
{HALF_INV_PHI, -0.5, -HALF_PHI, 0},
|
||||
{HALF_INV_PHI, 0, 0.5, -HALF_PHI},
|
||||
{HALF_INV_PHI, -HALF_PHI, 0, 0.5},
|
||||
{HALF_INV_PHI, 0.5, -HALF_PHI, 0},
|
||||
{HALF_INV_PHI, HALF_PHI, 0, -0.5},
|
||||
{HALF_INV_PHI, -0.5, HALF_PHI, 0},
|
||||
{HALF_INV_PHI, 0, -0.5, HALF_PHI},
|
||||
{HALF_INV_PHI, HALF_PHI, 0, 0.5},
|
||||
{HALF_INV_PHI, 0, 0.5, HALF_PHI},
|
||||
{HALF_INV_PHI, 0.5, HALF_PHI, 0},
|
||||
{0, 1, 0, 0},
|
||||
{0, HALF_PHI, -0.5, HALF_INV_PHI},
|
||||
{0, HALF_PHI, -0.5, -HALF_INV_PHI},
|
||||
{0, HALF_PHI, 0.5, HALF_INV_PHI},
|
||||
{0, HALF_PHI, 0.5, -HALF_INV_PHI},
|
||||
{0, 0.5, HALF_INV_PHI, -HALF_PHI},
|
||||
{0, 0.5, HALF_INV_PHI, HALF_PHI},
|
||||
{0, 0.5, -HALF_INV_PHI, -HALF_PHI},
|
||||
{0, 0.5, -HALF_INV_PHI, HALF_PHI},
|
||||
{0, HALF_INV_PHI, -HALF_PHI, 0.5},
|
||||
{0, HALF_INV_PHI, -HALF_PHI, -0.5},
|
||||
{0, HALF_INV_PHI, HALF_PHI, 0.5},
|
||||
{0, HALF_INV_PHI, HALF_PHI, -0.5},
|
||||
{0, 0, 1, 0},
|
||||
{0, 0, 0, 1} };
|
||||
|
||||
static void quat_rot(double* r, double* a, double* b)
|
||||
{
|
||||
b[0] = (r[0] * a[0] - r[1] * a[1] - r[2] * a[2] - r[3] * a[3]);
|
||||
b[1] = (r[0] * a[1] + r[1] * a[0] + r[2] * a[3] - r[3] * a[2]);
|
||||
b[2] = (r[0] * a[2] - r[1] * a[3] + r[2] * a[0] + r[3] * a[1]);
|
||||
b[3] = (r[0] * a[3] + r[1] * a[2] - r[2] * a[1] + r[3] * a[0]);
|
||||
}
|
||||
|
||||
static int rotate_quaternion_into_fundamental_zone(int num_generators, double (*generator)[4], double* q)
|
||||
{
|
||||
double max = 0.0;
|
||||
int i = 0, bi = -1;
|
||||
for (i=0;i<num_generators;i++)
|
||||
{
|
||||
double* g = generator[i];
|
||||
double t = fabs(q[0] * g[0] - q[1] * g[1] - q[2] * g[2] - q[3] * g[3]);
|
||||
if (t > max)
|
||||
{
|
||||
max = t;
|
||||
bi = i;
|
||||
}
|
||||
}
|
||||
|
||||
double f[4];
|
||||
quat_rot(q, generator[bi], f);
|
||||
memcpy(q, &f, 4 * sizeof(double));
|
||||
if (q[0] < 0)
|
||||
{
|
||||
q[0] = -q[0];
|
||||
q[1] = -q[1];
|
||||
q[2] = -q[2];
|
||||
q[3] = -q[3];
|
||||
}
|
||||
|
||||
return bi;
|
||||
}
|
||||
|
||||
int rotate_quaternion_into_cubic_fundamental_zone(double* q)
|
||||
{
|
||||
return rotate_quaternion_into_fundamental_zone(24, generator_cubic, q);
|
||||
}
|
||||
|
||||
int rotate_quaternion_into_diamond_cubic_fundamental_zone(double* q)
|
||||
{
|
||||
return rotate_quaternion_into_fundamental_zone(12, generator_diamond_cubic, q);
|
||||
}
|
||||
|
||||
int rotate_quaternion_into_icosahedral_fundamental_zone(double* q)
|
||||
{
|
||||
return rotate_quaternion_into_fundamental_zone(60, generator_icosahedral, q);
|
||||
}
|
||||
|
||||
int rotate_quaternion_into_hcp_fundamental_zone(double* q)
|
||||
{
|
||||
return rotate_quaternion_into_fundamental_zone(6, generator_hcp, q);
|
||||
}
|
||||
|
||||
int rotate_quaternion_into_diamond_hexagonal_fundamental_zone(double* q)
|
||||
{
|
||||
return rotate_quaternion_into_fundamental_zone(3, generator_diamond_hexagonal, q);
|
||||
}
|
||||
|
||||
double quat_dot(double* a, double* b)
|
||||
{
|
||||
return a[0] * b[0]
|
||||
+ a[1] * b[1]
|
||||
+ a[2] * b[2]
|
||||
+ a[3] * b[3];
|
||||
}
|
||||
|
||||
double quat_size(double* q)
|
||||
{
|
||||
return sqrt(quat_dot(q, q));
|
||||
}
|
||||
|
||||
void normalize_quaternion(double* q)
|
||||
{
|
||||
double size = quat_size(q);
|
||||
|
||||
q[0] /= size;
|
||||
q[1] /= size;
|
||||
q[2] /= size;
|
||||
q[3] /= size;
|
||||
}
|
||||
|
||||
void rotation_matrix_to_quaternion(double* u, double* q)
|
||||
{
|
||||
double r11 = u[0];
|
||||
double r12 = u[1];
|
||||
double r13 = u[2];
|
||||
double r21 = u[3];
|
||||
double r22 = u[4];
|
||||
double r23 = u[5];
|
||||
double r31 = u[6];
|
||||
double r32 = u[7];
|
||||
double r33 = u[8];
|
||||
|
||||
q[0] = (1.0 + r11 + r22 + r33) / 4.0;
|
||||
q[1] = (1.0 + r11 - r22 - r33) / 4.0;
|
||||
q[2] = (1.0 - r11 + r22 - r33) / 4.0;
|
||||
q[3] = (1.0 - r11 - r22 + r33) / 4.0;
|
||||
|
||||
q[0] = sqrt(MAX(0, q[0]));
|
||||
q[1] = sqrt(MAX(0, q[1]));
|
||||
q[2] = sqrt(MAX(0, q[2]));
|
||||
q[3] = sqrt(MAX(0, q[3]));
|
||||
|
||||
double m0 = MAX(q[0], q[1]);
|
||||
double m1 = MAX(q[2], q[3]);
|
||||
double max = MAX(m0, m1);
|
||||
|
||||
int i = 0;
|
||||
for (i=0;i<4;i++)
|
||||
if (q[i] == max)
|
||||
break;
|
||||
|
||||
if (i == 0)
|
||||
{
|
||||
q[1] *= SIGN(r32 - r23);
|
||||
q[2] *= SIGN(r13 - r31);
|
||||
q[3] *= SIGN(r21 - r12);
|
||||
}
|
||||
else if (i == 1)
|
||||
{
|
||||
q[0] *= SIGN(r32 - r23);
|
||||
q[2] *= SIGN(r21 + r12);
|
||||
q[3] *= SIGN(r13 + r31);
|
||||
}
|
||||
else if (i == 2)
|
||||
{
|
||||
q[0] *= SIGN(r13 - r31);
|
||||
q[1] *= SIGN(r21 + r12);
|
||||
q[3] *= SIGN(r32 + r23);
|
||||
}
|
||||
else if (i == 3)
|
||||
{
|
||||
q[0] *= SIGN(r21 - r12);
|
||||
q[1] *= SIGN(r31 + r13);
|
||||
q[2] *= SIGN(r32 + r23);
|
||||
}
|
||||
|
||||
normalize_quaternion(q);
|
||||
}
|
||||
|
||||
void quaternion_to_rotation_matrix(double* q, double* u)
|
||||
{
|
||||
double a = q[0];
|
||||
double b = q[1];
|
||||
double c = q[2];
|
||||
double d = q[3];
|
||||
|
||||
u[0] = a*a + b*b - c*c - d*d;
|
||||
u[1] = 2*b*c - 2*a*d;
|
||||
u[2] = 2*b*d + 2*a*c;
|
||||
|
||||
u[3] = 2*b*c + 2*a*d;
|
||||
u[4] = a*a - b*b + c*c - d*d;
|
||||
u[5] = 2*c*d - 2*a*b;
|
||||
|
||||
u[6] = 2*b*d - 2*a*c;
|
||||
u[7] = 2*c*d + 2*a*b;
|
||||
u[8] = a*a - b*b - c*c + d*d;
|
||||
}
|
||||
|
||||
double quat_quick_misorientation(double* q1, double* q2)
|
||||
{
|
||||
double t = quat_dot(q1, q2);
|
||||
t = MIN(1, MAX(-1, t));
|
||||
return 2 * t * t - 1;
|
||||
}
|
||||
|
||||
double quat_misorientation(double* q1, double* q2)
|
||||
{
|
||||
return acos(quat_quick_misorientation(q1, q2));
|
||||
}
|
||||
|
||||
|
||||
double quat_quick_disorientation_cubic(double* q0, double* q1)
|
||||
{
|
||||
double qrot[4];
|
||||
double qinv[4] = {q0[0], -q0[1], -q0[2], -q0[3]};
|
||||
quat_rot(qinv, q1, qrot);
|
||||
|
||||
rotate_quaternion_into_cubic_fundamental_zone(qrot);
|
||||
double t = qrot[0];
|
||||
t = MIN(1, MAX(-1, t));
|
||||
return 2 * t * t - 1;
|
||||
}
|
||||
|
||||
double quat_disorientation_cubic(double* q0, double* q1)
|
||||
{
|
||||
return acos(quat_quick_disorientation_cubic(q0, q1));
|
||||
}
|
||||
|
||||
double quat_quick_disorientation_diamond_cubic(double* q0, double* q1)
|
||||
{
|
||||
double qrot[4];
|
||||
double qinv[4] = {q0[0], -q0[1], -q0[2], -q0[3]};
|
||||
quat_rot(qinv, q1, qrot);
|
||||
|
||||
rotate_quaternion_into_diamond_cubic_fundamental_zone(qrot);
|
||||
double t = qrot[0];
|
||||
t = MIN(1, MAX(-1, t));
|
||||
return 2 * t * t - 1;
|
||||
}
|
||||
|
||||
double quat_disorientation_diamond_cubic(double* q0, double* q1)
|
||||
{
|
||||
return acos(quat_quick_disorientation_diamond_cubic(q0, q1));
|
||||
}
|
||||
|
||||
double quat_quick_disorientation_hcp(double* q0, double* q1)
|
||||
{
|
||||
double qrot[4];
|
||||
double qinv[4] = {q0[0], -q0[1], -q0[2], -q0[3]};
|
||||
quat_rot(qinv, q1, qrot);
|
||||
|
||||
rotate_quaternion_into_hcp_fundamental_zone(qrot);
|
||||
double t = qrot[0];
|
||||
t = MIN(1, MAX(-1, t));
|
||||
return 2 * t * t - 1;
|
||||
}
|
||||
|
||||
double quat_disorientation_hcp(double* q0, double* q1)
|
||||
{
|
||||
return acos(quat_quick_disorientation_hcp(q0, q1));
|
||||
}
|
||||
|
||||
double quat_quick_disorientation_diamond_hexagonal(double* q0, double* q1)
|
||||
{
|
||||
double qrot[4];
|
||||
double qinv[4] = {q0[0], -q0[1], -q0[2], -q0[3]};
|
||||
quat_rot(qinv, q1, qrot);
|
||||
|
||||
rotate_quaternion_into_diamond_hexagonal_fundamental_zone(qrot);
|
||||
double t = qrot[0];
|
||||
t = MIN(1, MAX(-1, t));
|
||||
return 2 * t * t - 1;
|
||||
}
|
||||
|
||||
double quat_disorientation_diamond_hexagonal(double* q0, double* q1)
|
||||
{
|
||||
return acos(quat_quick_disorientation_diamond_hexagonal(q0, q1));
|
||||
}
|
||||
|
||||
double quat_quick_disorientation_icosahedral(double* q0, double* q1)
|
||||
{
|
||||
double qrot[4];
|
||||
double qinv[4] = {q0[0], -q0[1], -q0[2], -q0[3]};
|
||||
quat_rot(qinv, q1, qrot);
|
||||
|
||||
rotate_quaternion_into_icosahedral_fundamental_zone(qrot);
|
||||
double t = qrot[0];
|
||||
t = MIN(1, MAX(-1, t));
|
||||
return 2 * t * t - 1;
|
||||
}
|
||||
|
||||
double quat_disorientation_icosahedral(double* q0, double* q1)
|
||||
{
|
||||
return acos(quat_quick_disorientation_icosahedral(q0, q1));
|
||||
}
|
||||
|
||||
32
src/PTM/ptm_quat.h
Normal file
32
src/PTM/ptm_quat.h
Normal file
@ -0,0 +1,32 @@
|
||||
#ifndef PTM_QUAT_H
|
||||
#define PTM_QUAT_H
|
||||
|
||||
int rotate_quaternion_into_cubic_fundamental_zone(double* q);
|
||||
int rotate_quaternion_into_diamond_cubic_fundamental_zone(double* q);
|
||||
int rotate_quaternion_into_icosahedral_fundamental_zone(double* q);
|
||||
int rotate_quaternion_into_hcp_fundamental_zone(double* q);
|
||||
int rotate_quaternion_into_diamond_hexagonal_fundamental_zone(double* q);
|
||||
|
||||
void normalize_quaternion(double* q);
|
||||
void quaternion_to_rotation_matrix(double* q, double* U);
|
||||
void rotation_matrix_to_quaternion(double* u, double* q);
|
||||
double quat_dot(double* a, double* b);
|
||||
double quat_quick_misorientation(double* q1, double* q2);
|
||||
double quat_misorientation(double* q1, double* q2);
|
||||
|
||||
double quat_quick_disorientation_cubic(double* q0, double* q1);
|
||||
double quat_disorientation_cubic(double* q0, double* q1);
|
||||
double quat_quick_disorientation_diamond_cubic(double* q0, double* q1);
|
||||
double quat_disorientation_diamond_cubic(double* q0, double* q1);
|
||||
double quat_quick_disorientation_hcp(double* q0, double* q1);
|
||||
double quat_disorientation_hcp(double* q0, double* q1);
|
||||
double quat_quick_disorientation_diamond_hexagonal(double* q0, double* q1);
|
||||
double quat_disorientation_diamond_hexagonal(double* q0, double* q1);
|
||||
double quat_quick_disorientation_icosahedral(double* q0, double* q1);
|
||||
double quat_disorientation_icosahedral(double* q0, double* q1);
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
294
src/PTM/ptm_structure_matcher.cpp
Normal file
294
src/PTM/ptm_structure_matcher.cpp
Normal file
@ -0,0 +1,294 @@
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <string.h>
|
||||
#include <cmath>
|
||||
#include <cfloat>
|
||||
#include <cassert>
|
||||
#include <algorithm>
|
||||
#include "ptm_convex_hull_incremental.h"
|
||||
#include "ptm_canonical_coloured.h"
|
||||
#include "ptm_graph_data.h"
|
||||
#include "ptm_graph_tools.h"
|
||||
#include "ptm_normalize_vertices.h"
|
||||
#include "ptm_polar.h"
|
||||
#include "ptm_structure_matcher.h"
|
||||
#include "ptm_constants.h"
|
||||
|
||||
|
||||
static double calc_rmsd(int num_points, const double (*ideal_points)[3], double (*normalized)[3], int8_t* mapping,
|
||||
double G1, double G2, double E0, double* q, double* p_scale)
|
||||
{
|
||||
double A0[9];
|
||||
InnerProduct(A0, num_points, ideal_points, normalized, mapping);
|
||||
|
||||
double nrmsdsq, rot[9];
|
||||
FastCalcRMSDAndRotation(A0, E0, &nrmsdsq, q, rot);
|
||||
|
||||
double k0 = 0;
|
||||
for (int i=0;i<num_points;i++)
|
||||
{
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
double v = 0.0;
|
||||
for (int k=0;k<3;k++)
|
||||
v += rot[j*3+k] * ideal_points[i][k];
|
||||
|
||||
k0 += v * normalized[mapping[i]][j];
|
||||
}
|
||||
}
|
||||
|
||||
double scale = k0 / G2;
|
||||
*p_scale = scale;
|
||||
return sqrt(fabs(G1 - scale*k0) / num_points);
|
||||
}
|
||||
|
||||
static void check_graphs( const refdata_t* s,
|
||||
uint64_t hash,
|
||||
int8_t* canonical_labelling,
|
||||
double (*normalized)[3],
|
||||
result_t* res)
|
||||
{
|
||||
int num_points = s->num_nbrs + 1;
|
||||
const double (*ideal_points)[3] = s->points;
|
||||
int8_t inverse_labelling[PTM_MAX_POINTS];
|
||||
int8_t mapping[PTM_MAX_POINTS];
|
||||
|
||||
for (int i=0; i<num_points; i++)
|
||||
inverse_labelling[ canonical_labelling[i] ] = i;
|
||||
|
||||
double G1 = 0, G2 = 0;
|
||||
for (int i=0;i<num_points;i++)
|
||||
{
|
||||
double x1 = ideal_points[i][0];
|
||||
double y1 = ideal_points[i][1];
|
||||
double z1 = ideal_points[i][2];
|
||||
|
||||
double x2 = normalized[i][0];
|
||||
double y2 = normalized[i][1];
|
||||
double z2 = normalized[i][2];
|
||||
|
||||
G1 += x1 * x1 + y1 * y1 + z1 * z1;
|
||||
G2 += x2 * x2 + y2 * y2 + z2 * z2;
|
||||
}
|
||||
double E0 = (G1 + G2) / 2;
|
||||
|
||||
for (int i = 0;i<s->num_graphs;i++)
|
||||
{
|
||||
if (hash != s->graphs[i].hash)
|
||||
continue;
|
||||
|
||||
graph_t* gref = &s->graphs[i];
|
||||
for (int j = 0;j<gref->num_automorphisms;j++)
|
||||
{
|
||||
for (int k=0;k<num_points;k++)
|
||||
mapping[automorphisms[gref->automorphism_index + j][k]] = inverse_labelling[ gref->canonical_labelling[k] ];
|
||||
|
||||
double q[4], scale = 0;
|
||||
double rmsd = calc_rmsd(num_points, ideal_points, normalized, mapping, G1, G2, E0, q, &scale);
|
||||
if (rmsd < res->rmsd)
|
||||
{
|
||||
res->rmsd = rmsd;
|
||||
res->scale = scale;
|
||||
res->ref_struct = s;
|
||||
memcpy(res->q, q, 4 * sizeof(double));
|
||||
memcpy(res->mapping, mapping, sizeof(int8_t) * num_points);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int match_general(const refdata_t* s, double (*ch_points)[3], double (*points)[3], convexhull_t* ch, result_t* res)
|
||||
{
|
||||
int8_t degree[PTM_MAX_NBRS];
|
||||
int8_t facets[PTM_MAX_FACETS][3];
|
||||
|
||||
int ret = get_convex_hull(s->num_nbrs + 1, (const double (*)[3])ch_points, ch, facets);
|
||||
ch->ok = ret >= 0;
|
||||
if (ret != 0)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
if (ch->num_facets != s->num_facets)
|
||||
return PTM_NO_ERROR; //incorrect number of facets in convex hull
|
||||
|
||||
int max_degree = graph_degree(s->num_facets, facets, s->num_nbrs, degree);
|
||||
if (max_degree > s->max_degree)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
if (s->type == PTM_MATCH_SC)
|
||||
for (int i = 0;i<s->num_nbrs;i++)
|
||||
if (degree[i] != 4)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
double normalized[PTM_MAX_POINTS][3];
|
||||
subtract_barycentre(s->num_nbrs + 1, points, normalized);
|
||||
|
||||
int8_t code[2 * PTM_MAX_EDGES];
|
||||
int8_t colours[PTM_MAX_POINTS] = {0};
|
||||
int8_t canonical_labelling[PTM_MAX_POINTS];
|
||||
uint64_t hash = 0;
|
||||
ret = canonical_form_coloured(s->num_facets, facets, s->num_nbrs, degree, colours, canonical_labelling, &code[0], &hash);
|
||||
if (ret != PTM_NO_ERROR)
|
||||
return ret;
|
||||
|
||||
check_graphs(s, hash, canonical_labelling, normalized, res);
|
||||
return PTM_NO_ERROR;
|
||||
}
|
||||
|
||||
int match_fcc_hcp_ico(double (*ch_points)[3], double (*points)[3], int32_t flags, convexhull_t* ch, result_t* res)
|
||||
{
|
||||
int num_nbrs = structure_fcc.num_nbrs;
|
||||
int num_facets = structure_fcc.num_facets;
|
||||
int max_degree = structure_fcc.max_degree;
|
||||
|
||||
int8_t degree[PTM_MAX_NBRS];
|
||||
int8_t facets[PTM_MAX_FACETS][3];
|
||||
|
||||
int ret = get_convex_hull(num_nbrs + 1, (const double (*)[3])ch_points, ch, facets);
|
||||
ch->ok = ret >= 0;
|
||||
if (ret != 0)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
if (ch->num_facets != num_facets)
|
||||
return PTM_NO_ERROR; //incorrect number of facets in convex hull
|
||||
|
||||
int _max_degree = graph_degree(num_facets, facets, num_nbrs, degree);
|
||||
if (_max_degree > max_degree)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
double normalized[PTM_MAX_POINTS][3];
|
||||
subtract_barycentre(num_nbrs + 1, points, normalized);
|
||||
|
||||
int8_t code[2 * PTM_MAX_EDGES];
|
||||
int8_t colours[PTM_MAX_POINTS] = {0};
|
||||
int8_t canonical_labelling[PTM_MAX_POINTS];
|
||||
uint64_t hash = 0;
|
||||
ret = canonical_form_coloured(num_facets, facets, num_nbrs, degree, colours, canonical_labelling, &code[0], &hash);
|
||||
if (ret != PTM_NO_ERROR)
|
||||
return ret;
|
||||
|
||||
if (flags & PTM_CHECK_FCC) check_graphs(&structure_fcc, hash, canonical_labelling, normalized, res);
|
||||
if (flags & PTM_CHECK_HCP) check_graphs(&structure_hcp, hash, canonical_labelling, normalized, res);
|
||||
if (flags & PTM_CHECK_ICO) check_graphs(&structure_ico, hash, canonical_labelling, normalized, res);
|
||||
return PTM_NO_ERROR;
|
||||
}
|
||||
|
||||
int match_dcub_dhex(double (*ch_points)[3], double (*points)[3], int32_t flags, convexhull_t* ch, result_t* res)
|
||||
{
|
||||
int num_nbrs = structure_dcub.num_nbrs;
|
||||
int num_facets = structure_fcc.num_facets;
|
||||
int max_degree = structure_dcub.max_degree;
|
||||
|
||||
|
||||
int8_t facets[PTM_MAX_FACETS][3];
|
||||
int ret = get_convex_hull(num_nbrs + 1, (const double (*)[3])ch_points, ch, facets);
|
||||
ch->ok = ret >= 0;
|
||||
if (ret != 0)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
//check for facets with multiple inner atoms
|
||||
bool inverted[4] = {false, false, false, false};
|
||||
for (int i=0;i<ch->num_facets;i++)
|
||||
{
|
||||
int n = 0;
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
if (facets[i][j] <= 3)
|
||||
{
|
||||
inverted[facets[i][j]] = true;
|
||||
n++;
|
||||
}
|
||||
}
|
||||
if (n > 1)
|
||||
return PTM_NO_ERROR;
|
||||
}
|
||||
|
||||
int num_inverted = 0;
|
||||
for (int i=0;i<4;i++)
|
||||
num_inverted += inverted[i] ? 1 : 0;
|
||||
|
||||
if (ch->num_facets != num_facets + 2 * num_inverted)
|
||||
return PTM_NO_ERROR; //incorrect number of facets in convex hull
|
||||
|
||||
int8_t degree[PTM_MAX_NBRS];
|
||||
int _max_degree = graph_degree(num_facets, facets, num_nbrs, degree);
|
||||
if (_max_degree > max_degree)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
int num_found = 0;
|
||||
int8_t toadd[4][3];
|
||||
for (int i=0;i<ch->num_facets;i++)
|
||||
{
|
||||
int a = facets[i][0];
|
||||
int b = facets[i][1];
|
||||
int c = facets[i][2];
|
||||
if (a <= 3 || b <= 3 || c <= 3)
|
||||
continue;
|
||||
|
||||
int i0 = (a - 4) / 3;
|
||||
int i1 = (b - 4) / 3;
|
||||
int i2 = (c - 4) / 3;
|
||||
|
||||
if (i0 == i1 && i0 == i2)
|
||||
{
|
||||
if (num_found + num_inverted >= 4)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
toadd[num_found][0] = a;
|
||||
toadd[num_found][1] = b;
|
||||
toadd[num_found][2] = c;
|
||||
num_found++;
|
||||
|
||||
memcpy(&facets[i], &facets[ch->num_facets - 1], 3 * sizeof(int8_t));
|
||||
ch->num_facets--;
|
||||
i--;
|
||||
}
|
||||
}
|
||||
|
||||
if (num_found + num_inverted != 4)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
for (int i=0;i<num_found;i++)
|
||||
{
|
||||
int a = toadd[i][0];
|
||||
int b = toadd[i][1];
|
||||
int c = toadd[i][2];
|
||||
|
||||
int i0 = (a - 4) / 3;
|
||||
|
||||
facets[ch->num_facets][0] = i0;
|
||||
facets[ch->num_facets][1] = b;
|
||||
facets[ch->num_facets][2] = c;
|
||||
ch->num_facets++;
|
||||
|
||||
facets[ch->num_facets][0] = a;
|
||||
facets[ch->num_facets][1] = i0;
|
||||
facets[ch->num_facets][2] = c;
|
||||
ch->num_facets++;
|
||||
|
||||
facets[ch->num_facets][0] = a;
|
||||
facets[ch->num_facets][1] = b;
|
||||
facets[ch->num_facets][2] = i0;
|
||||
ch->num_facets++;
|
||||
}
|
||||
|
||||
_max_degree = graph_degree(ch->num_facets, facets, num_nbrs, degree);
|
||||
if (_max_degree > max_degree)
|
||||
return PTM_NO_ERROR;
|
||||
|
||||
double normalized[PTM_MAX_POINTS][3];
|
||||
subtract_barycentre(num_nbrs + 1, points, normalized);
|
||||
|
||||
int8_t code[2 * PTM_MAX_EDGES];
|
||||
int8_t colours[PTM_MAX_POINTS] = {1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
|
||||
int8_t canonical_labelling[PTM_MAX_POINTS];
|
||||
uint64_t hash = 0;
|
||||
ret = canonical_form_coloured(ch->num_facets, facets, num_nbrs, degree, colours, canonical_labelling, &code[0], &hash);
|
||||
if (ret != PTM_NO_ERROR)
|
||||
return ret;
|
||||
|
||||
if (flags & PTM_CHECK_DCUB) check_graphs(&structure_dcub, hash, canonical_labelling, normalized, res);
|
||||
if (flags & PTM_CHECK_DHEX) check_graphs(&structure_dhex, hash, canonical_labelling, normalized, res);
|
||||
|
||||
return PTM_NO_ERROR;
|
||||
}
|
||||
|
||||
21
src/PTM/ptm_structure_matcher.h
Normal file
21
src/PTM/ptm_structure_matcher.h
Normal file
@ -0,0 +1,21 @@
|
||||
#ifndef PTM_STRUCTURE_MATCHER_H
|
||||
#define PTM_STRUCTURE_MATCHER_H
|
||||
|
||||
#include "ptm_initialize_data.h"
|
||||
#include "ptm_constants.h"
|
||||
|
||||
typedef struct
|
||||
{
|
||||
double rmsd;
|
||||
double scale;
|
||||
double q[4]; //rotation in quaternion form (rigid body transformation)
|
||||
int8_t mapping[PTM_MAX_POINTS];
|
||||
const refdata_t* ref_struct;
|
||||
} result_t;
|
||||
|
||||
int match_general(const refdata_t* s, double (*ch_points)[3], double (*points)[3], convexhull_t* ch, result_t* res);
|
||||
int match_fcc_hcp_ico(double (*ch_points)[3], double (*points)[3], int32_t flags, convexhull_t* ch, result_t* res);
|
||||
int match_dcub_dhex(double (*ch_points)[3], double (*points)[3], int32_t flags, convexhull_t* ch, result_t* res);
|
||||
|
||||
#endif
|
||||
|
||||
1368
src/PTM/ptm_voronoi_cell.cpp
Normal file
1368
src/PTM/ptm_voronoi_cell.cpp
Normal file
File diff suppressed because it is too large
Load Diff
324
src/PTM/ptm_voronoi_cell.h
Normal file
324
src/PTM/ptm_voronoi_cell.h
Normal file
@ -0,0 +1,324 @@
|
||||
// Voro++, a 3D cell-based Voronoi library
|
||||
//
|
||||
// Author : Chris H. Rycroft (LBL / UC Berkeley)
|
||||
// Email : chr@alum.mit.edu
|
||||
// Date : August 30th 2011
|
||||
//
|
||||
// Modified by PM Larsen for use in Polyhedral Template Matching
|
||||
|
||||
/** \file cell.hh
|
||||
* \brief Header file for the voronoicell and related classes. */
|
||||
|
||||
#ifndef PTM_VOROPP_CELL_HH
|
||||
#define PTM_VOROPP_CELL_HH
|
||||
|
||||
#include <vector>
|
||||
#include <cstdio>
|
||||
|
||||
#include "ptm_voronoi_config.h"
|
||||
|
||||
namespace voro {
|
||||
|
||||
/** \brief A class representing a single Voronoi cell.
|
||||
*
|
||||
* This class represents a single Voronoi cell, as a collection of vertices
|
||||
* that are connected by edges. The class contains routines for initializing
|
||||
* the Voronoi cell to be simple shapes such as a box, tetrahedron, or octahedron.
|
||||
* It the contains routines for recomputing the cell based on cutting it
|
||||
* by a plane, which forms the key routine for the Voronoi cell computation.
|
||||
* It contains numerous routine for computing statistics about the Voronoi cell,
|
||||
* and it can output the cell in several formats.
|
||||
*
|
||||
* This class is not intended for direct use, but forms the base of the
|
||||
* voronoicell and voronoicell_neighbor classes, which extend it based on
|
||||
* whether neighboring particle ID information needs to be tracked. */
|
||||
class voronoicell_base {
|
||||
public:
|
||||
/** This holds the current size of the arrays ed and nu, which
|
||||
* hold the vertex information. If more vertices are created
|
||||
* than can fit in this array, then it is dynamically extended
|
||||
* using the add_memory_vertices routine. */
|
||||
int current_vertices;
|
||||
/** This holds the current maximum allowed order of a vertex,
|
||||
* which sets the size of the mem, mep, and mec arrays. If a
|
||||
* vertex is created with more vertices than this, the arrays
|
||||
* are dynamically extended using the add_memory_vorder routine.
|
||||
*/
|
||||
int current_vertex_order;
|
||||
/** This sets the size of the main delete stack. */
|
||||
int current_delete_size;
|
||||
/** This sets the size of the auxiliary delete stack. */
|
||||
int current_delete2_size;
|
||||
/** This sets the total number of vertices in the current cell.
|
||||
*/
|
||||
int p;
|
||||
/** This is the index of particular point in the cell, which is
|
||||
* used to start the tracing routines for plane intersection
|
||||
* and cutting. These routines will work starting from any
|
||||
* point, but it's often most efficient to start from the last
|
||||
* point considered, since in many cases, the cell construction
|
||||
* algorithm may consider many planes with similar vectors
|
||||
* concurrently. */
|
||||
int up;
|
||||
/** This is a two dimensional array that holds information
|
||||
* about the edge connections of the vertices that make up the
|
||||
* cell. The two dimensional array is not allocated in the
|
||||
* usual method. To account for the fact the different vertices
|
||||
* have different orders, and thus require different amounts of
|
||||
* storage, the elements of ed[i] point to one-dimensional
|
||||
* arrays in the mep[] array of different sizes.
|
||||
*
|
||||
* More specifically, if vertex i has order m, then ed[i]
|
||||
* points to a one-dimensional array in mep[m] that has 2*m+1
|
||||
* entries. The first m elements hold the neighboring edges, so
|
||||
* that the jth edge of vertex i is held in ed[i][j]. The next
|
||||
* m elements hold a table of relations which is redundant but
|
||||
* helps speed up the computation. It satisfies the relation
|
||||
* ed[ed[i][j]][ed[i][m+j]]=i. The final entry holds a back
|
||||
* pointer, so that ed[i+2*m]=i. The back pointers are used
|
||||
* when rearranging the memory. */
|
||||
int **ed;
|
||||
/** This array holds the order of the vertices in the Voronoi
|
||||
* cell. This array is dynamically allocated, with its current
|
||||
* size held by current_vertices. */
|
||||
int *nu;
|
||||
/** This in an array with size 3*current_vertices for holding
|
||||
* the positions of the vertices. */
|
||||
double *pts;
|
||||
voronoicell_base();
|
||||
virtual ~voronoicell_base();
|
||||
void init_base(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax);
|
||||
void init_octahedron_base(double l);
|
||||
void init_tetrahedron_base(double x0,double y0,double z0,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3);
|
||||
void translate(double x,double y,double z);
|
||||
double volume();
|
||||
double max_radius_squared();
|
||||
double total_edge_distance();
|
||||
double surface_area();
|
||||
void centroid(double &cx,double &cy,double &cz);
|
||||
int number_of_faces();
|
||||
int number_of_edges();
|
||||
void vertex_orders(std::vector<int> &v);
|
||||
void vertices(std::vector<double> &v);
|
||||
void vertices(double x,double y,double z,std::vector<double> &v);
|
||||
void face_areas(std::vector<double> &v);
|
||||
void face_orders(std::vector<int> &v);
|
||||
void face_freq_table(std::vector<int> &v);
|
||||
void face_vertices(std::vector<int> &v);
|
||||
void face_perimeters(std::vector<double> &v);
|
||||
void normals(std::vector<double> &v);
|
||||
template<class vc_class>
|
||||
bool nplane(vc_class &vc,double x,double y,double z,double rsq,int p_id);
|
||||
bool plane_intersects(double x,double y,double z,double rsq);
|
||||
bool plane_intersects_guess(double x,double y,double z,double rsq);
|
||||
void construct_relations();
|
||||
void check_relations();
|
||||
void check_duplicates();
|
||||
/** Returns a list of IDs of neighboring particles
|
||||
* corresponding to each face.
|
||||
* \param[out] v a reference to a vector in which to return the
|
||||
* results. If no neighbor information is
|
||||
* available, a blank vector is returned. */
|
||||
virtual void neighbors(std::vector<int> &v) {v.clear();}
|
||||
/** This a virtual function that is overridden by a routine to
|
||||
* print the neighboring particle IDs for a given vertex. By
|
||||
* default, when no neighbor information is available, the
|
||||
* routine does nothing.
|
||||
* \param[in] i the vertex to consider. */
|
||||
/** This is a simple inline function for picking out the index
|
||||
* of the next edge counterclockwise at the current vertex.
|
||||
* \param[in] a the index of an edge of the current vertex.
|
||||
* \param[in] p the number of the vertex.
|
||||
* \return 0 if a=nu[p]-1, or a+1 otherwise. */
|
||||
inline int cycle_up(int a,int p) {return a==nu[p]-1?0:a+1;}
|
||||
/** This is a simple inline function for picking out the index
|
||||
* of the next edge clockwise from the current vertex.
|
||||
* \param[in] a the index of an edge of the current vertex.
|
||||
* \param[in] p the number of the vertex.
|
||||
* \return nu[p]-1 if a=0, or a-1 otherwise. */
|
||||
inline int cycle_down(int a,int p) {return a==0?nu[p]-1:a-1;}
|
||||
protected:
|
||||
/** This a one dimensional array that holds the current sizes
|
||||
* of the memory allocations for them mep array.*/
|
||||
int *mem;
|
||||
/** This is a one dimensional array that holds the current
|
||||
* number of vertices of order p that are stored in the mep[p]
|
||||
* array. */
|
||||
int *mec;
|
||||
/** This is a two dimensional array for holding the information
|
||||
* about the edges of the Voronoi cell. mep[p] is a
|
||||
* one-dimensional array for holding the edge information about
|
||||
* all vertices of order p, with each vertex holding 2*p+1
|
||||
* integers of information. The total number of vertices held
|
||||
* on mep[p] is stored in mem[p]. If the space runs out, the
|
||||
* code allocates more using the add_memory() routine. */
|
||||
int **mep;
|
||||
inline void reset_edges();
|
||||
template<class vc_class>
|
||||
void check_memory_for_copy(vc_class &vc,voronoicell_base* vb);
|
||||
void copy(voronoicell_base* vb);
|
||||
private:
|
||||
/** This is the delete stack, used to store the vertices which
|
||||
* are going to be deleted during the plane cutting procedure.
|
||||
*/
|
||||
int *ds,*stacke;
|
||||
/** This is the auxiliary delete stack, which has size set by
|
||||
* current_delete2_size. */
|
||||
int *ds2,*stacke2;
|
||||
/** This stores the current memory allocation for the marginal
|
||||
* cases. */
|
||||
int current_marginal;
|
||||
/** This stores the total number of marginal points which are
|
||||
* currently in the buffer. */
|
||||
int n_marg;
|
||||
/** This array contains a list of the marginal points, and also
|
||||
* the outcomes of the marginal tests. */
|
||||
int *marg;
|
||||
/** The x coordinate of the normal vector to the test plane. */
|
||||
double px;
|
||||
/** The y coordinate of the normal vector to the test plane. */
|
||||
double py;
|
||||
/** The z coordinate of the normal vector to the test plane. */
|
||||
double pz;
|
||||
/** The magnitude of the normal vector to the test plane. */
|
||||
double prsq;
|
||||
template<class vc_class>
|
||||
void add_memory(vc_class &vc,int i,int *stackp2);
|
||||
template<class vc_class>
|
||||
void add_memory_vertices(vc_class &vc);
|
||||
template<class vc_class>
|
||||
void add_memory_vorder(vc_class &vc);
|
||||
void add_memory_ds(int *&stackp);
|
||||
void add_memory_ds2(int *&stackp2);
|
||||
template<class vc_class>
|
||||
inline bool collapse_order1(vc_class &vc);
|
||||
template<class vc_class>
|
||||
inline bool collapse_order2(vc_class &vc);
|
||||
template<class vc_class>
|
||||
inline bool delete_connection(vc_class &vc,int j,int k,bool hand);
|
||||
template<class vc_class>
|
||||
inline bool search_for_outside_edge(vc_class &vc,int &up);
|
||||
template<class vc_class>
|
||||
inline void add_to_stack(vc_class &vc,int lp,int *&stackp2);
|
||||
inline bool plane_intersects_track(double x,double y,double z,double rs,double g);
|
||||
inline void normals_search(std::vector<double> &v,int i,int j,int k);
|
||||
inline bool search_edge(int l,int &m,int &k);
|
||||
inline int m_test(int n,double &ans);
|
||||
int check_marginal(int n,double &ans);
|
||||
friend class voronoicell;
|
||||
friend class voronoicell_neighbor;
|
||||
};
|
||||
|
||||
/** \brief Extension of the voronoicell_base class to represent a Voronoi cell
|
||||
* with neighbor information.
|
||||
*
|
||||
* This class is an extension of the voronoicell_base class, in cases when the
|
||||
* IDs of neighboring particles associated with each face of the Voronoi cell.
|
||||
* It contains additional data structures mne and ne for storing this
|
||||
* information. */
|
||||
class voronoicell_neighbor : public voronoicell_base {
|
||||
public:
|
||||
using voronoicell_base::nplane;
|
||||
/** This two dimensional array holds the neighbor information
|
||||
* associated with each vertex. mne[p] is a one dimensional
|
||||
* array which holds all of the neighbor information for
|
||||
* vertices of order p. */
|
||||
int **mne;
|
||||
/** This is a two dimensional array that holds the neighbor
|
||||
* information associated with each vertex. ne[i] points to a
|
||||
* one-dimensional array in mne[nu[i]]. ne[i][j] holds the
|
||||
* neighbor information associated with the jth edge of vertex
|
||||
* i. It is set to the ID number of the plane that made the
|
||||
* face that is clockwise from the jth edge. */
|
||||
int **ne;
|
||||
voronoicell_neighbor();
|
||||
~voronoicell_neighbor();
|
||||
void operator=(voronoicell_neighbor &c);
|
||||
/** Cuts the Voronoi cell by a particle whose center is at a
|
||||
* separation of (x,y,z) from the cell center. The value of rsq
|
||||
* should be initially set to \f$x^2+y^2+z^2\f$.
|
||||
* \param[in] (x,y,z) the normal vector to the plane.
|
||||
* \param[in] rsq the distance along this vector of the plane.
|
||||
* \param[in] p_id the plane ID (for neighbor tracking only).
|
||||
* \return False if the plane cut deleted the cell entirely,
|
||||
* true otherwise. */
|
||||
inline bool nplane(double x,double y,double z,double rsq,int p_id) {
|
||||
return nplane(*this,x,y,z,rsq,p_id);
|
||||
}
|
||||
/** This routine calculates the modulus squared of the vector
|
||||
* before passing it to the main nplane() routine with full
|
||||
* arguments.
|
||||
* \param[in] (x,y,z) the vector to cut the cell by.
|
||||
* \param[in] p_id the plane ID (for neighbor tracking only).
|
||||
* \return False if the plane cut deleted the cell entirely,
|
||||
* true otherwise. */
|
||||
inline bool nplane(double x,double y,double z,int p_id) {
|
||||
double rsq=x*x+y*y+z*z;
|
||||
return nplane(*this,x,y,z,rsq,p_id);
|
||||
}
|
||||
/** This version of the plane routine just makes up the plane
|
||||
* ID to be zero. It will only be referenced if neighbor
|
||||
* tracking is enabled.
|
||||
* \param[in] (x,y,z) the vector to cut the cell by.
|
||||
* \param[in] rsq the modulus squared of the vector.
|
||||
* \return False if the plane cut deleted the cell entirely,
|
||||
* true otherwise. */
|
||||
inline bool plane(double x,double y,double z,double rsq) {
|
||||
return nplane(*this,x,y,z,rsq,0);
|
||||
}
|
||||
/** Cuts a Voronoi cell using the influence of a particle at
|
||||
* (x,y,z), first calculating the modulus squared of this
|
||||
* vector before passing it to the main nplane() routine. Zero
|
||||
* is supplied as the plane ID, which will be ignored unless
|
||||
* neighbor tracking is enabled.
|
||||
* \param[in] (x,y,z) the vector to cut the cell by.
|
||||
* \return False if the plane cut deleted the cell entirely,
|
||||
* true otherwise. */
|
||||
inline bool plane(double x,double y,double z) {
|
||||
double rsq=x*x+y*y+z*z;
|
||||
return nplane(*this,x,y,z,rsq,0);
|
||||
}
|
||||
void init(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax);
|
||||
void check_facets();
|
||||
virtual void neighbors(std::vector<int> &v);
|
||||
|
||||
private:
|
||||
int *paux1;
|
||||
int *paux2;
|
||||
inline void n_allocate(int i,int m) {mne[i]=new int[m*i];}
|
||||
inline void n_add_memory_vertices(int i) {
|
||||
int **pp=new int*[i];
|
||||
for(int j=0;j<current_vertices;j++) pp[j]=ne[j];
|
||||
delete [] ne;ne=pp;
|
||||
}
|
||||
inline void n_add_memory_vorder(int i) {
|
||||
int **p2=new int*[i];
|
||||
for(int j=0;j<current_vertex_order;j++) p2[j]=mne[j];
|
||||
delete [] mne;mne=p2;
|
||||
}
|
||||
inline void n_set_pointer(int p,int n) {
|
||||
ne[p]=mne[n]+n*mec[n];
|
||||
}
|
||||
inline void n_copy(int a,int b,int c,int d) {ne[a][b]=ne[c][d];}
|
||||
inline void n_set(int a,int b,int c) {ne[a][b]=c;}
|
||||
inline void n_set_aux1(int k) {paux1=mne[k]+k*mec[k];}
|
||||
inline void n_copy_aux1(int a,int b) {paux1[b]=ne[a][b];}
|
||||
inline void n_copy_aux1_shift(int a,int b) {paux1[b]=ne[a][b+1];}
|
||||
inline void n_set_aux2_copy(int a,int b) {
|
||||
paux2=mne[b]+b*mec[b];
|
||||
for(int i=0;i<b;i++) ne[a][i]=paux2[i];
|
||||
}
|
||||
inline void n_copy_pointer(int a,int b) {ne[a]=ne[b];}
|
||||
inline void n_set_to_aux1(int j) {ne[j]=paux1;}
|
||||
inline void n_set_to_aux2(int j) {ne[j]=paux2;}
|
||||
inline void n_allocate_aux1(int i) {paux1=new int[i*mem[i]];}
|
||||
inline void n_switch_to_aux1(int i) {delete [] mne[i];mne[i]=paux1;}
|
||||
inline void n_copy_to_aux1(int i,int m) {paux1[m]=mne[i][m];}
|
||||
inline void n_set_to_aux1_offset(int k,int m) {ne[k]=paux1+m;}
|
||||
friend class voronoicell_base;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
129
src/PTM/ptm_voronoi_config.h
Normal file
129
src/PTM/ptm_voronoi_config.h
Normal file
@ -0,0 +1,129 @@
|
||||
// Voro++, a 3D cell-based Voronoi library
|
||||
//
|
||||
// Author : Chris H. Rycroft (LBL / UC Berkeley)
|
||||
// Email : chr@alum.mit.edu
|
||||
// Date : August 30th 2011
|
||||
//
|
||||
// Modified by PM Larsen for use in Polyhedral Template Matching
|
||||
|
||||
/** \file config.hh
|
||||
* \brief Master configuration file for setting various compile-time options. */
|
||||
|
||||
#ifndef PTM_VOROPP_CONFIG_HH
|
||||
#define PTM_VOROPP_CONFIG_HH
|
||||
|
||||
namespace voro {
|
||||
|
||||
// These constants set the initial memory allocation for the Voronoi cell
|
||||
/** The initial memory allocation for the number of vertices. */
|
||||
const int init_vertices=256;
|
||||
/** The initial memory allocation for the maximum vertex order. */
|
||||
const int init_vertex_order=64;
|
||||
/** The initial memory allocation for the number of regular vertices of order
|
||||
* 3. */
|
||||
const int init_3_vertices=256;
|
||||
/** The initial memory allocation for the number of vertices of higher order.
|
||||
*/
|
||||
const int init_n_vertices=8;
|
||||
/** The initial buffer size for marginal cases used by the suretest class. */
|
||||
const int init_marginal=64;
|
||||
/** The initial size for the delete stack. */
|
||||
const int init_delete_size=256;
|
||||
/** The initial size for the auxiliary delete stack. */
|
||||
const int init_delete2_size=256;
|
||||
/** The initial size for the wall pointer array. */
|
||||
const int init_wall_size=32;
|
||||
/** The default initial size for the ordering class. */
|
||||
const int init_ordering_size=4096;
|
||||
/** The initial size of the pre_container chunk index. */
|
||||
const int init_chunk_size=256;
|
||||
|
||||
// If the initial memory is too small, the program dynamically allocates more.
|
||||
// However, if the limits below are reached, then the program bails out.
|
||||
/** The maximum memory allocation for the number of vertices. */
|
||||
const int max_vertices=16777216;
|
||||
/** The maximum memory allocation for the maximum vertex order. */
|
||||
const int max_vertex_order=2048;
|
||||
/** The maximum memory allocation for the any particular order of vertex. */
|
||||
const int max_n_vertices=16777216;
|
||||
/** The maximum buffer size for marginal cases used by the suretest class. */
|
||||
const int max_marginal=16777216;
|
||||
/** The maximum size for the delete stack. */
|
||||
const int max_delete_size=16777216;
|
||||
/** The maximum size for the auxiliary delete stack. */
|
||||
const int max_delete2_size=16777216;
|
||||
/** The maximum amount of particle memory allocated for a single region. */
|
||||
const int max_particle_memory=16777216;
|
||||
/** The maximum size for the wall pointer array. */
|
||||
const int max_wall_size=2048;
|
||||
/** The maximum size for the ordering class. */
|
||||
const int max_ordering_size=67108864;
|
||||
/** The maximum size for the pre_container chunk index. */
|
||||
const int max_chunk_size=65536;
|
||||
|
||||
/** The chunk size in the pre_container classes. */
|
||||
const int pre_container_chunk_size=1024;
|
||||
|
||||
#ifndef VOROPP_VERBOSE
|
||||
/** Voro++ can print a number of different status and debugging messages to
|
||||
* notify the user of special behavior, and this macro sets the amount which
|
||||
* are displayed. At level 0, no messages are printed. At level 1, messages
|
||||
* about unusual cases during cell construction are printed, such as when the
|
||||
* plane routine bails out due to floating point problems. At level 2, general
|
||||
* messages about memory expansion are printed. At level 3, technical details
|
||||
* about memory management are printed. */
|
||||
#define VOROPP_VERBOSE 0
|
||||
#endif
|
||||
|
||||
/** If a point is within this distance of a cutting plane, then the code
|
||||
* assumes that point exactly lies on the plane. */
|
||||
const double tolerance=1e-11;
|
||||
|
||||
/** If a point is within this distance of a cutting plane, then the code stores
|
||||
* whether this point is inside, outside, or exactly on the cutting plane in
|
||||
* the marginal cases buffer, to prevent the test giving a different result on
|
||||
* a subsequent evaluation due to floating point rounding errors. */
|
||||
const double tolerance2=2e-11;
|
||||
|
||||
/** The square of the tolerance, used when deciding whether some squared
|
||||
* quantities are large enough to be used. */
|
||||
const double tolerance_sq=tolerance*tolerance;
|
||||
|
||||
/** A large number that is used in the computation. */
|
||||
const double large_number=1e30;
|
||||
|
||||
/** A radius to use as a placeholder when no other information is available. */
|
||||
const double default_radius=0.5;
|
||||
|
||||
/** The maximum number of shells of periodic images to test over. */
|
||||
const int max_unit_voro_shells=10;
|
||||
|
||||
/** A guess for the optimal number of particles per block, used to set up the
|
||||
* container grid. */
|
||||
const double optimal_particles=5.6;
|
||||
|
||||
/** If this is set to 1, then the code reports any instances of particles being
|
||||
* put outside of the container geometry. */
|
||||
#define VOROPP_REPORT_OUT_OF_BOUNDS 0
|
||||
|
||||
/** Voro++ returns this status code if there is a file-related error, such as
|
||||
* not being able to open file. */
|
||||
#define VOROPP_FILE_ERROR 1
|
||||
|
||||
/** Voro++ returns this status code if there is a memory allocation error, if
|
||||
* one of the safe memory limits is exceeded. */
|
||||
#define VOROPP_MEMORY_ERROR 2
|
||||
|
||||
/** Voro++ returns this status code if there is any type of internal error, if
|
||||
* it detects that representation of the Voronoi cell is inconsistent. This
|
||||
* status code will generally indicate a bug, and the developer should be
|
||||
* contacted. */
|
||||
#define VOROPP_INTERNAL_ERROR 3
|
||||
|
||||
/** Voro++ returns this status code if it could not interpret the command line
|
||||
* arguments passed to the command line utility. */
|
||||
#define VOROPP_CMD_LINE_ERROR 4
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user