tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixed folder structure

Browse Source
This commit is contained in:
pmla
2018-09-20 00:04:07 -04:00
parent c705e8d0e6
commit 37201beda5
63 changed files with 1 additions and 7293 deletions
Download Patch File Download Diff File Expand all files Collapse all files

307
src/PTM/compute_ptm_atom.cpp
View File

@ -1,307 +0,0 @@
/* ----------------------------------------------------------------------
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
View File

@ -1,48 +0,0 @@
/* -*- 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

174
src/PTM/ptm_constants.h
View File

@ -1,174 +0,0 @@
#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

27
src/PTM/ptm_functions.h
View File

@ -1,27 +0,0 @@
#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

101
src/USER-PTM/alloy_types.cpp
View File

@ -1,101 +0,0 @@
#include <algorithm>
#include "ptm_constants.h"
#include "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/USER-PTM/alloy_types.h
View File

@ -1,9 +0,0 @@
#ifndef ALLOY_TYPES_H
#define ALLOY_TYPES_H
#include "initialize_data.h"
int32_t find_alloy_type(const refdata_t* ref, int8_t* mapping, int32_t* numbers);
#endif

167
src/USER-PTM/canonical_coloured.cpp
View File

@ -1,167 +0,0 @@
#include <string.h>
#include <climits>
#include <algorithm>
#include "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/USER-PTM/canonical_coloured.h
View File

@ -1,9 +0,0 @@
#ifndef CANONICAL_COLOURED_H
#define 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

1368
src/USER-PTM/cell.cpp
View File

File diff suppressed because it is too large Load Diff

324
src/USER-PTM/cell.h
View File

@ -1,324 +0,0 @@
// 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 VOROPP_CELL_HH
#define VOROPP_CELL_HH
#include <vector>
#include <cstdio>
#include "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/USER-PTM/config.h
View File

@ -1,129 +0,0 @@
// 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 VOROPP_CONFIG_HH
#define 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

363
src/USER-PTM/convex_hull_incremental.cpp
View File

@ -1,363 +0,0 @@
#include <cmath>
#include <cfloat>
#include <string.h>
#include <cassert>
#include <algorithm>
#include "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/USER-PTM/convex_hull_incremental.h
View File

@ -1,27 +0,0 @@
#ifndef CONVEX_HULL_INCREMENTAL_H
#define 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/USER-PTM/deformation_gradient.cpp
View File

@ -1,37 +0,0 @@
#include "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/USER-PTM/deformation_gradient.h
View File

@ -1,142 +0,0 @@
#ifndef DEFORMATION_GRADIENT_H
#define 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

180
src/USER-PTM/fundamental_mappings.h
View File

@ -1,180 +0,0 @@
#ifndef FUNDAMENTAL_MAPPINGS_H
#define 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/USER-PTM/graph_data.cpp
View File

File diff suppressed because it is too large Load Diff

37
src/USER-PTM/graph_data.h
View File

@ -1,37 +0,0 @@
#ifndef GRAPH_DATA_H
#define 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/USER-PTM/graph_tools.cpp
View File

@ -1,52 +0,0 @@
#include <string.h>
#include <algorithm>
#include "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/USER-PTM/graph_tools.h
View File

@ -1,11 +0,0 @@
#ifndef GRAPH_TOOLS_H
#define 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/USER-PTM/index_ptm.cpp
View File

@ -1,218 +0,0 @@
#include <cstdio>
#include <cstdlib>
#include <string.h>
#include <cmath>
#include <cfloat>
#include <cassert>
#include <algorithm>
#include "convex_hull_incremental.h"
#include "graph_data.h"
#include "deformation_gradient.h"
#include "alloy_types.h"
#include "neighbour_ordering.h"
#include "normalize_vertices.h"
#include "quat.h"
#include "polar.h"
#include "initialize_data.h"
#include "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/USER-PTM/initialize_data.cpp
View File

@ -1,71 +0,0 @@
#include <cstdio>
#include <cstdlib>
#include <string.h>
#include <cmath>
#include <cfloat>
#include <cassert>
#include <algorithm>
#include "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/USER-PTM/initialize_data.h
View File

@ -1,61 +0,0 @@
#ifndef INITIALIZE_DATA_H
#define INITIALIZE_DATA_H
#include "graph_data.h"
#include "graph_tools.h"
#include "deformation_gradient.h"
#include "fundamental_mappings.h"
#include "neighbour_ordering.h"
#include "canonical_coloured.h"
#include "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/USER-PTM/neighbour_ordering.cpp
View File

@ -1,203 +0,0 @@
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <cassert>
#include <algorithm>
#include "ptm_constants.h"
#include "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/USER-PTM/neighbour_ordering.h
View File

@ -1,13 +0,0 @@
#ifndef NEIGHBOUR_ORDERING_H
#define 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/USER-PTM/normalize_vertices.cpp
View File

@ -1,55 +0,0 @@
#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/USER-PTM/normalize_vertices.h
View File

@ -1,8 +0,0 @@
#ifndef NORMALIZE_VERTICES_H
#define 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/USER-PTM/polar.cpp
View File

@ -1,337 +0,0 @@
/*******************************************************************************
* -/_|:|_|_\-
*
* 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 "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/USER-PTM/polar.h
View File

@ -1,12 +0,0 @@
#ifndef POLAR_H
#define 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

0
src/PTM/ptm_alloy_types.cpp → src/USER-PTM/ptm_alloy_types.cpp
View File

0
src/PTM/ptm_alloy_types.h → src/USER-PTM/ptm_alloy_types.h
View File

0
src/PTM/ptm_canonical_coloured.cpp → src/USER-PTM/ptm_canonical_coloured.cpp
View File

0
src/PTM/ptm_canonical_coloured.h → src/USER-PTM/ptm_canonical_coloured.h
View File

0
src/PTM/ptm_convex_hull_incremental.cpp → src/USER-PTM/ptm_convex_hull_incremental.cpp
View File

0
src/PTM/ptm_convex_hull_incremental.h → src/USER-PTM/ptm_convex_hull_incremental.h
View File

0
src/PTM/ptm_deformation_gradient.cpp → src/USER-PTM/ptm_deformation_gradient.cpp
View File

0
src/PTM/ptm_deformation_gradient.h → src/USER-PTM/ptm_deformation_gradient.h
View File

2
src/USER-PTM/ptm_functions.h
View File

@ -3,7 +3,7 @@
#include <stdint.h>
#include <stdbool.h>
#include "initialize_data.h"
#include "ptm_initialize_data.h"
#include "ptm_constants.h"

0
src/PTM/ptm_fundamental_mappings.h → src/USER-PTM/ptm_fundamental_mappings.h
View File

0
src/PTM/ptm_graph_data.cpp → src/USER-PTM/ptm_graph_data.cpp
View File

0
src/PTM/ptm_graph_data.h → src/USER-PTM/ptm_graph_data.h
View File

0
src/PTM/ptm_graph_tools.cpp → src/USER-PTM/ptm_graph_tools.cpp
View File

0
src/PTM/ptm_graph_tools.h → src/USER-PTM/ptm_graph_tools.h
View File

0
src/PTM/ptm_index.cpp → src/USER-PTM/ptm_index.cpp
View File

0
src/PTM/ptm_initialize_data.cpp → src/USER-PTM/ptm_initialize_data.cpp
View File

0
src/PTM/ptm_initialize_data.h → src/USER-PTM/ptm_initialize_data.h
View File

0
src/PTM/ptm_neighbour_ordering.cpp → src/USER-PTM/ptm_neighbour_ordering.cpp
View File

0
src/PTM/ptm_neighbour_ordering.h → src/USER-PTM/ptm_neighbour_ordering.h
View File

0
src/PTM/ptm_normalize_vertices.cpp → src/USER-PTM/ptm_normalize_vertices.cpp
View File

0
src/PTM/ptm_normalize_vertices.h → src/USER-PTM/ptm_normalize_vertices.h
View File

0
src/PTM/ptm_polar.cpp → src/USER-PTM/ptm_polar.cpp
View File

0
src/PTM/ptm_polar.h → src/USER-PTM/ptm_polar.h
View File

0
src/PTM/ptm_quat.cpp → src/USER-PTM/ptm_quat.cpp
View File

0
src/PTM/ptm_quat.h → src/USER-PTM/ptm_quat.h
View File

0
src/PTM/ptm_structure_matcher.cpp → src/USER-PTM/ptm_structure_matcher.cpp
View File

0
src/PTM/ptm_structure_matcher.h → src/USER-PTM/ptm_structure_matcher.h
View File

0
src/PTM/ptm_voronoi_cell.cpp → src/USER-PTM/ptm_voronoi_cell.cpp
View File

0
src/PTM/ptm_voronoi_cell.h → src/USER-PTM/ptm_voronoi_cell.h
View File

0
src/PTM/ptm_voronoi_config.h → src/USER-PTM/ptm_voronoi_config.h
View File

396
src/USER-PTM/quat.cpp
View File

@ -1,396 +0,0 @@
#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/USER-PTM/quat.h
View File

@ -1,32 +0,0 @@
#ifndef QUAT_H
#define 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/USER-PTM/structure_matcher.cpp
View File

@ -1,294 +0,0 @@
#include <cstdio>
#include <cstdlib>
#include <string.h>
#include <cmath>
#include <cfloat>
#include <cassert>
#include <algorithm>
#include "convex_hull_incremental.h"
#include "canonical_coloured.h"
#include "graph_data.h"
#include "graph_tools.h"
#include "normalize_vertices.h"
#include "polar.h"
#include "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/USER-PTM/structure_matcher.h
View File

@ -1,21 +0,0 @@
#ifndef STRUCTURE_MATCHER_H
#define STRUCTURE_MATCHER_H
#include "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