Merge pull request #1339 from pmla/ptm-update

sc (simple cubic) = 5

dcub (diamond cubic) = 6

dhex (diamond hexagonal) = 7

other = 8 :ul

graphene = 8 :ul

The value of the PTM structure will be 0 for atoms not in the specified

The value of the PTM structure will be 0 for unknown types and -1 for atoms not in the specified

compute group.  The choice of structures to search for can be specified using the "structures"

argument, which is a hyphen-separated list of structure keywords.

Two convenient pre-set options are provided:

default: fcc-hcp-bcc-ico

all: fcc-hcp-bcc-ico-sc-dcub-dhex :ul

all: fcc-hcp-bcc-ico-sc-dcub-dhex-graphene :ul

The 'default' setting detects the same structures as the Common Neighbor Analysis method.

The 'all' setting searches for all structure types.  A small performance penalty is

incurred for the diamond structures, so it is not recommended to use this option if

it is known that the simulation does not contain diamond structures.

The 'all' setting searches for all structure types.  A performance penalty is

incurred for the diamond and graphene structures, so it is not recommended to use this option if

it is known that the simulation does not contain these structures.

PTM identifies structures using two steps.  First, a graph isomorphism test is used

qw

qx

qy

qw :ul

qz :ul

The type is a number from 0 to 8.  The rmsd is a positive real number.

The type is a number from -1 to 8.  The rmsd is a positive real number.

The interatomic distance is computed from the scale factor in the RMSD equation.

The (qw,qx,qy,qz) parameters represent the orientation of the local structure

in quaternion form.  The reference coordinates for each template (from which the

# LAMMPS Input File for Grain Boundaries

# Mark Tschopp, Dec2009

# Modified by PM Larsen to demonstrate the use of PTM in LAMMPS

# This file will generate a single Sigma5(310) STGB and run PTM

# ---------- Initialize Simulation ---------------------

clear

units metal

dimension 3

boundary p p p

atom_style atomic

# ---------- Create Sigma5(310) Structure ---------

lattice fcc 4.05

region whole block 0.000000 12.807225 -64.0361225 64.0361225 0.000000 4.050000 units box

create_box 2 whole

region upper block INF INF 0.000000 64.0361225 INF INF units box

lattice fcc 4.05 orient x  0  3  1 orient y  0 -1  3 orient z  1  0  0

create_atoms 1 region upper

region lower block INF INF -64.0361225 0.000000 INF INF units box

lattice fcc 4.05 orient x  0  3 -1 orient y  0  1  3 orient z  1  0  0

create_atoms 2 region lower

group upper type 1

group lower type 2

mass 1 1.0

mass 2 1.0

# ---------- Define Interatomic Potential ---------------------

pair_style	lj/cut 2.5

pair_coeff	* * 1 1

pair_coeff	1 1 1 1.1 2.8

# ---------- Displace atoms and delete overlapping atoms ---------------------

displace_atoms upper move 0 0 0 units lattice

delete_atoms overlap 0.35 lower upper

# ---------- Define PTM settings (default structures, RMSD threshold of 0.1) -------------------

compute ptm all ptm/atom default 0.1

# ---------- Dump data into Data file -------------

reset_timestep 0

dump 		1 all cfg 10000 dump.ptm_example_*.cfg mass type xs ys zs c_ptm[1] c_ptm[2] c_ptm[3] c_ptm[4] c_ptm[5] c_ptm[6] c_ptm[7]

dump_modify     1 element Al Al

run 1

print "All done"

http://arxiv.org/abs/1603.05143

Contact email: pmla@mit.edu

#include <cmath>

#include <cstdlib>

#include <cstring>

#include <vector>

#include "atom.h"

#include "comm.h"

#include "ptm_functions.h"

#define MAX_NEIGHBORS 30

#define NUM_COLUMNS 7

#define UNKNOWN 0

#define OTHER 8

#define PTM_LAMMPS_UNKNOWN -1

#define PTM_LAMMPS_OTHER 0

using namespace LAMMPS_NS;

  char *ptr = structures;

  const char *strings[] = {"fcc",  "hcp",  "bcc", "ico",    "sc",

                           "dcub", "dhex", "all", "default"};

                           "dcub", "dhex", "graphene", "all", "default"};

  int num_strings = sizeof(strings) / sizeof(const char*);

  int32_t flags[] = {

      PTM_CHECK_FCC,

      PTM_CHECK_HCP,

      PTM_CHECK_SC,

      PTM_CHECK_DCUB,

      PTM_CHECK_DHEX,

      PTM_CHECK_GRAPHENE,

      PTM_CHECK_ALL,

      PTM_CHECK_FCC | PTM_CHECK_HCP | PTM_CHECK_BCC | PTM_CHECK_ICO};

  while (*ptr != '\0') {

    bool found = false;

    for (int i = 0; i < 9; i++) {

    for (int i = 0; i < num_strings; i++) {

      int len = strlen(strings[i]);

      if (strncmp(ptr, strings[i], len) == 0) {

        input_flags |= flags[i];

/* ---------------------------------------------------------------------- */

void ComputePTMAtom::init_list(int id, NeighList *ptr) { list = ptr; }

void ComputePTMAtom::init_list(int /* id */, NeighList *ptr) { list = ptr; }

/* ---------------------------------------------------------------------- */

typedef struct

{

  double **x;

  int *numneigh;

  int **firstneigh;

  int *ilist;

  int nlocal;

} ptmnbrdata_t;

typedef struct {

  int index;

  double d;

  return a.d < b.d;

}

static int get_neighbors(double *pos, int jnum, int *jlist, double **x,

                         double (*nbr)[3]) {

static int get_neighbours(void* vdata, size_t central_index, size_t atom_index, int num, size_t* nbr_indices, int32_t* numbers, double (*nbr_pos)[3])

{

  ptmnbrdata_t* data = (ptmnbrdata_t*)vdata;

  double **x = data->x;

  double *pos = x[atom_index];

  int *jlist = NULL;

  int jnum = 0;

  if (atom_index < data->nlocal) {

    jlist = data->firstneigh[atom_index];

    jnum = data->numneigh[atom_index];

  }

  else {

    jlist = data->firstneigh[central_index];

    jnum = data->numneigh[central_index];

  }

  ptmnbr_t *nbr_order = new ptmnbr_t[jnum];

  std::vector<ptmnbr_t> nbr_order;

  for (int jj = 0; jj < jnum; jj++) {

    int j = jlist[jj];

    j &= NEIGHMASK;

    if (j == atom_index)

      continue;

    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;

    ptmnbr_t nbr = {j, rsq};

    nbr_order.push_back(nbr);

  }

  std::sort(nbr_order, nbr_order + jnum, &sorthelper_compare);

  int num_nbrs = std::min(MAX_NEIGHBORS, jnum);

  std::sort(nbr_order.begin(), nbr_order.end(), &sorthelper_compare);

  int num_nbrs = std::min(num - 1, (int)nbr_order.size());

  nbr[0][0] = nbr[0][1] = nbr[0][2] = 0;

  nbr_pos[0][0] = nbr_pos[0][1] = nbr_pos[0][2] = 0;

  nbr_indices[0] = atom_index;

  numbers[0] = 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];

    nbr_pos[jj + 1][0] = x[j][0] - pos[0];

    nbr_pos[jj + 1][1] = x[j][1] - pos[1];

    nbr_pos[jj + 1][2] = x[j][2] - pos[2];

    nbr_indices[jj + 1] = j;

    numbers[jj + 1] = 0;

  }

  delete[] nbr_order;

  return num_nbrs;

  return num_nbrs + 1;

}

void ComputePTMAtom::compute_peratom() {

  double **x = atom->x;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  ptmnbrdata_t nbrlist = {x, numneigh, firstneigh, ilist, atom->nlocal};

  for (int ii = 0; ii < inum; ii++) {

    int i = ilist[ii];

    output[i][0] = UNKNOWN;

    output[i][0] = PTM_LAMMPS_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);

    double scale, rmsd, interatomic_distance;

    double q[4];

    bool standard_orientations = false;

    ptm_index(local_handle, i, get_neighbours, (void*)&nbrlist,

              input_flags, standard_orientations,

              &type, &alloy_type, &scale, &rmsd, q,

              NULL, NULL, NULL, NULL, &interatomic_distance, NULL, NULL);

    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;

    if (type == PTM_MATCH_NONE) {

      type = PTM_LAMMPS_OTHER;

      rmsd = INFINITY;

    }

    output[i][0] = type;

    output[i][1] = rmsd;

/*Copyright (c) 2016 PM Larsen

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

*/

#include <algorithm>

#include "ptm_constants.h"

#include "ptm_initialize_data.h"

#define NUM_ALLOY_TYPES 3

static uint32_t typedata[NUM_ALLOY_TYPES][3] = {

	{PTM_MATCH_FCC, PTM_ALLOY_L10,    0x000001fe},

	{PTM_MATCH_FCC, PTM_ALLOY_L12_CU, 0x0000001e},

        {PTM_MATCH_FCC, PTM_ALLOY_L10,    0x00000db6},

        {PTM_MATCH_FCC, PTM_ALLOY_L12_CU, 0x00000492},

        {PTM_MATCH_FCC, PTM_ALLOY_L12_AU, 0x00001ffe},

};

int32_t find_alloy_type(const refdata_t* ref, int8_t* mapping, int32_t* numbers)

{

        for (int i=0;i<ref->num_nbrs+1;i++)

                if (numbers[i] == -1)

                        return PTM_ALLOY_NONE;

        if (test_pure(ref->num_nbrs, numbers))

                return PTM_ALLOY_PURE;

                if (test_shell_structure(ref, mapping, numbers, 4))

                        return PTM_ALLOY_SIC;

        if (ref->type == PTM_MATCH_GRAPHENE)

                if (test_shell_structure(ref, mapping, numbers, 3))

                        return PTM_ALLOY_BN;

        return PTM_ALLOY_NONE;

}