diff --git a/src/MISC/fix_orient_bcc.cpp b/src/MISC/fix_orient_bcc.cpp new file mode 100644 index 0000000000..88db82bd82 --- /dev/null +++ b/src/MISC/fix_orient_bcc.cpp @@ -0,0 +1,606 @@ +/* ---------------------------------------------------------------------- + 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 authors: Koenraad Janssens and David Olmsted (SNL) + Modification for bcc provided by: Tegar Wicaksono (UBC) + For a tutorial, please see "Order parameters of crystals in LAMMPS" + (https://dx.doi.org/10.6084/m9.figshare.1488628.v1 +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "fix_orient_bcc.h" +#include "atom.h" +#include "update.h" +#include "respa.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "comm.h" +#include "output.h" +#include "force.h" +#include "math_const.h" +#include "citeme.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathConst; + +#define BIG 1000000000 + +static const char cite_fix_orient_bcc[] = + "fix orient/bcc command:\n\n" + "@Article{Wicaksono16,\n" + " author = {A. T. Wicaksono, C. W. Sinclair, M. Militzer},\n" + " title = {An atomistic study of the correlation between the migration of planar and curved grain boundaries},\n" + " journal = {Computational Materials Science},\n" + " year = 2016,\n" + " volume = 117,\n" + " pages = {397--405}\n" + "}\n\n"; + +/* ---------------------------------------------------------------------- */ + +FixOrientBCC::FixOrientBCC(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (lmp->citeme) lmp->citeme->add(cite_fix_orient_bcc); + + MPI_Comm_rank(world,&me); + + if (narg != 11) error->all(FLERR,"Illegal fix orient/bcc command"); + + scalar_flag = 1; + global_freq = 1; + extscalar = 1; + + peratom_flag = 1; + size_peratom_cols = 2; + peratom_freq = 1; + + nstats = force->inumeric(FLERR,arg[3]); + direction_of_motion = force->inumeric(FLERR,arg[4]); + a = force->numeric(FLERR,arg[5]); + Vxi = force->numeric(FLERR,arg[6]); + uxif_low = force->numeric(FLERR,arg[7]); + uxif_high = force->numeric(FLERR,arg[8]); + + if (direction_of_motion == 0) { + int n = strlen(arg[9]) + 1; + chifilename = new char[n]; + strcpy(chifilename,arg[9]); + n = strlen(arg[10]) + 1; + xifilename = new char[n]; + strcpy(xifilename,arg[10]); + } else if (direction_of_motion == 1) { + int n = strlen(arg[9]) + 1; + xifilename = new char[n]; + strcpy(xifilename,arg[9]); + n = strlen(arg[10]) + 1; + chifilename = new char[n]; + strcpy(chifilename,arg[10]); + } else error->all(FLERR,"Illegal fix orient/bcc command"); + + // initializations + + half_bcc_nn = 4; + use_xismooth = false; + double xicutoff = 1.57; + xicutoffsq = xicutoff * xicutoff; + cutsq = 0.75 * a*a*xicutoffsq; + nmax = 0; + + // read xi and chi reference orientations from files + + if (me == 0) { + char line[IMGMAX]; + char *result; + int count; + + FILE *infile = fopen(xifilename,"r"); + if (infile == NULL) error->one(FLERR,"Fix orient/bcc file open failed"); + for (int i = 0; i < 4; i++) { + result = fgets(line,IMGMAX,infile); + if (!result) error->one(FLERR,"Fix orient/bcc file read failed"); + count = sscanf(line,"%lg %lg %lg",&Rxi[i][0],&Rxi[i][1],&Rxi[i][2]); + if (count != 3) error->one(FLERR,"Fix orient/bcc file read failed"); + } + fclose(infile); + + infile = fopen(chifilename,"r"); + if (infile == NULL) error->one(FLERR,"Fix orient/bcc file open failed"); + for (int i = 0; i < 4; i++) { + result = fgets(line,IMGMAX,infile); + if (!result) error->one(FLERR,"Fix orient/bcc file read failed"); + count = sscanf(line,"%lg %lg %lg",&Rchi[i][0],&Rchi[i][1],&Rchi[i][2]); + if (count != 3) error->one(FLERR,"Fix orient/bcc file read failed"); + } + fclose(infile); + } + + MPI_Bcast(&Rxi[0][0],18,MPI_DOUBLE,0,world); + MPI_Bcast(&Rchi[0][0],18,MPI_DOUBLE,0,world); + + // make copy of the reference vectors + + for (int i = 0; i < 4; i++) + for (int j = 0; j < 3; j++) { + half_xi_chi_vec[0][i][j] = Rxi[i][j]; + half_xi_chi_vec[1][i][j] = Rchi[i][j]; + } + + // compute xiid,xi0,xi1 for all 8 neighbors + // xi is the favored crystal + // want order parameter when actual is Rchi + + double xi_sq,dxi[3],rchi[3]; + + xiid = 0.0; + for (int i = 0; i < 4; i++) { + rchi[0] = Rchi[i][0]; + rchi[1] = Rchi[i][1]; + rchi[2] = Rchi[i][2]; + find_best_ref(rchi,0,xi_sq,dxi); + xiid += sqrt(xi_sq); + for (int j = 0; j < 3; j++) rchi[j] = -rchi[j]; + find_best_ref(rchi,0,xi_sq,dxi); + xiid += sqrt(xi_sq); + } + + xiid /= 8.0; + xi0 = uxif_low * xiid; + xi1 = uxif_high * xiid; + + // set comm size needed by this Fix + // NOTE: doesn't seem that use_xismooth is ever true + + if (use_xismooth) comm_forward = 62; + else comm_forward = 50; + + added_energy = 0.0; + + nmax = atom->nmax; + nbr = (Nbr *) memory->smalloc(nmax*sizeof(Nbr),"orient/bcc:nbr"); + memory->create(order,nmax,2,"orient/bcc:order"); + array_atom = order; + + // zero the array since a variable may access it before first run + + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) order[i][0] = order[i][1] = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +FixOrientBCC::~FixOrientBCC() +{ + delete [] xifilename; + delete [] chifilename; + memory->sfree(nbr); + memory->destroy(order); +} + +/* ---------------------------------------------------------------------- */ + +int FixOrientBCC::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= THERMO_ENERGY; + mask |= POST_FORCE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientBCC::init() +{ + if (strstr(update->integrate_style,"respa")) + nlevels_respa = ((Respa *) update->integrate)->nlevels; + + // need a full neighbor list, built whenever re-neighboring occurs + + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->fix = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientBCC::init_list(int id, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientBCC::setup(int vflag) +{ + if (strstr(update->integrate_style,"verlet")) + post_force(vflag); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(vflag,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientBCC::post_force(int vflag) +{ + int i,j,k,ii,jj,inum,jnum,m,n,nn,nsort; + tagint id_self; + int *ilist,*jlist,*numneigh,**firstneigh; + double edelta,omega; + double dx,dy,dz,rsq,xismooth,xi_sq,duxi,duxi_other; + double dxi[3]; + double *dxiptr; + bool found_myself; + + // set local ptrs + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + tagint *tag = atom->tag; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // insure nbr and order data structures are adequate size + + if (nall > nmax) { + nmax = nall; + memory->destroy(nbr); + memory->destroy(order); + nbr = (Nbr *) memory->smalloc(nmax*sizeof(Nbr),"orient/bcc:nbr"); + memory->create(order,nmax,2,"orient/bcc:order"); + array_atom = order; + } + + // loop over owned atoms and build Nbr data structure of neighbors + // use full neighbor list + + added_energy = 0.0; + int count = 0; + int mincount = BIG; + int maxcount = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + if (jnum < mincount) mincount = jnum; + if (jnum > maxcount) { + if (maxcount) delete [] sort; + sort = new Sort[jnum]; + maxcount = jnum; + } + + // loop over all neighbors of atom i + // for those within cutsq, build sort data structure + // store local id, rsq, delta vector, xismooth (if included) + + nsort = 0; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + count++; + + dx = x[i][0] - x[j][0]; + dy = x[i][1] - x[j][1]; + dz = x[i][2] - x[j][2]; + rsq = dx*dx + dy*dy + dz*dz; + + if (rsq < cutsq) { + sort[nsort].id = j; + sort[nsort].rsq = rsq; + sort[nsort].delta[0] = dx; + sort[nsort].delta[1] = dy; + sort[nsort].delta[2] = dz; + if (use_xismooth) { + xismooth = (xicutoffsq - 4.0*rsq/(3.0*a*a)) / (xicutoffsq - 1.0); + sort[nsort].xismooth = 1.0 - fabs(1.0-xismooth); + } + nsort++; + } + } + + // sort neighbors by rsq distance + // no need to sort if nsort <= 8 + + if (nsort > 8) qsort(sort,nsort,sizeof(Sort),compare); + + // copy up to 8 nearest neighbors into nbr data structure + // operate on delta vector via find_best_ref() to compute dxi + + n = MIN(8,nsort); + nbr[i].n = n; + if (n == 0) continue; + + double xi_total = 0.0; + for (j = 0; j < n; j++) { + find_best_ref(sort[j].delta,0,xi_sq,dxi); + xi_total += sqrt(xi_sq); + nbr[i].id[j] = sort[j].id; + nbr[i].dxi[j][0] = dxi[0]/n; + nbr[i].dxi[j][1] = dxi[1]/n; + nbr[i].dxi[j][2] = dxi[2]/n; + if (use_xismooth) nbr[i].xismooth[j] = sort[j].xismooth; + } + xi_total /= n; + order[i][0] = xi_total; + + // compute potential derivative to xi + + if (xi_total < xi0) { + nbr[i].duxi = 0.0; + edelta = 0.0; + order[i][1] = 0.0; + } else if (xi_total > xi1) { + nbr[i].duxi = 0.0; + edelta = Vxi; + order[i][1] = 1.0; + } else { + omega = MY_PI2*(xi_total-xi0) / (xi1-xi0); + nbr[i].duxi = MY_PI*Vxi*sin(2.0*omega) / (2.0*(xi1-xi0)); + edelta = Vxi*(1 - cos(2.0*omega)) / 2.0; + order[i][1] = omega / MY_PI2; + } + added_energy += edelta; + } + + if (maxcount) delete [] sort; + + // communicate to acquire nbr data for ghost atoms + + comm->forward_comm_fix(this); + + // compute grain boundary force on each owned atom + // skip atoms not in group + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + n = nbr[i].n; + duxi = nbr[i].duxi; + + for (j = 0; j < n; j++) { + dxiptr = &nbr[i].dxi[j][0]; + if (use_xismooth) { + xismooth = nbr[i].xismooth[j]; + f[i][0] += duxi * dxiptr[0] * xismooth; + f[i][1] += duxi * dxiptr[1] * xismooth; + f[i][2] += duxi * dxiptr[2] * xismooth; + } else { + f[i][0] += duxi * dxiptr[0]; + f[i][1] += duxi * dxiptr[1]; + f[i][2] += duxi * dxiptr[2]; + } + + // m = local index of neighbor + // id_self = ID for atom I in atom M's neighbor list + // if M is local atom, id_self will be local ID of atom I + // if M is ghost atom, id_self will be global ID of atom I + + m = nbr[i].id[j]; + if (m < nlocal) id_self = i; + else id_self = tag[i]; + found_myself = false; + nn = nbr[m].n; + + for (k = 0; k < nn; k++) { + if (id_self == nbr[m].id[k]) { + if (found_myself) error->one(FLERR,"Fix orient/bcc found self twice"); + found_myself = true; + duxi_other = nbr[m].duxi; + dxiptr = &nbr[m].dxi[k][0]; + if (use_xismooth) { + xismooth = nbr[m].xismooth[k]; + f[i][0] -= duxi_other * dxiptr[0] * xismooth; + f[i][1] -= duxi_other * dxiptr[1] * xismooth; + f[i][2] -= duxi_other * dxiptr[2] * xismooth; + } else { + f[i][0] -= duxi_other * dxiptr[0]; + f[i][1] -= duxi_other * dxiptr[1]; + f[i][2] -= duxi_other * dxiptr[2]; + } + } + } + } + } + + // print statistics every nstats timesteps + + if (nstats && update->ntimestep % nstats == 0) { + int total; + MPI_Allreduce(&count,&total,1,MPI_INT,MPI_SUM,world); + double ave = static_cast(total)/atom->natoms; + + int min,max; + MPI_Allreduce(&mincount,&min,1,MPI_INT,MPI_MIN,world); + MPI_Allreduce(&maxcount,&max,1,MPI_INT,MPI_MAX,world); + + if (me == 0) { + if (screen) fprintf(screen, + "orient step " BIGINT_FORMAT ": " BIGINT_FORMAT + " atoms have %d neighbors\n", + update->ntimestep,atom->natoms,total); + if (logfile) fprintf(logfile, + "orient step " BIGINT_FORMAT ": " BIGINT_FORMAT + " atoms have %d neighbors\n", + update->ntimestep,atom->natoms,total); + if (screen) + fprintf(screen," neighs: min = %d, max = %d, ave = %g\n", + min,max,ave); + if (logfile) + fprintf(logfile," neighs: min = %d, max = %d, ave = %g\n", + min,max,ave); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientBCC::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +double FixOrientBCC::compute_scalar() +{ + double added_energy_total; + MPI_Allreduce(&added_energy,&added_energy_total,1,MPI_DOUBLE,MPI_SUM,world); + return added_energy_total; +} + +/* ---------------------------------------------------------------------- */ + +int FixOrientBCC::pack_forward_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,k,num; + tagint id; + + tagint *tag = atom->tag; + int nlocal = atom->nlocal; + int m = 0; + + for (i = 0; i < n; i++) { + k = list[i]; + num = nbr[k].n; + buf[m++] = num; + buf[m++] = nbr[k].duxi; + + for (j = 0; j < num; j++) { + if (use_xismooth) buf[m++] = nbr[k].xismooth[j]; + buf[m++] = nbr[k].dxi[j][0]; + buf[m++] = nbr[k].dxi[j][1]; + buf[m++] = nbr[k].dxi[j][2]; + + // id stored in buf needs to be global ID + // if k is a local atom, it stores local IDs, so convert to global + // if k is a ghost atom (already comm'd), its IDs are already global + + id = nbr[k].id[j]; + if (k < nlocal) id = tag[id]; + buf[m++] = id; + } + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientBCC::unpack_forward_comm(int n, int first, double *buf) +{ + int i,j,num; + int last = first + n; + int m = 0; + + for (i = first; i < last; i++) { + nbr[i].n = num = static_cast (buf[m++]); + nbr[i].duxi = buf[m++]; + + for (j = 0; j < num; j++) { + if (use_xismooth) nbr[i].xismooth[j] = buf[m++]; + nbr[i].dxi[j][0] = buf[m++]; + nbr[i].dxi[j][1] = buf[m++]; + nbr[i].dxi[j][2] = buf[m++]; + nbr[i].id[j] = static_cast (buf[m++]); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientBCC::find_best_ref(double *displs, int which_crystal, + double &xi_sq, double *dxi) +{ + int i; + double dot,tmp; + + double best_dot = -1.0; // best is biggest (smallest angle) + int best_i = -1; + int best_sign = 0; + + for (i = 0; i < half_bcc_nn; i++) { + dot = displs[0] * half_xi_chi_vec[which_crystal][i][0] + + displs[1] * half_xi_chi_vec[which_crystal][i][1] + + displs[2] * half_xi_chi_vec[which_crystal][i][2]; + if (fabs(dot) > best_dot) { + best_dot = fabs(dot); + best_i = i; + if (dot < 0.0) best_sign = -1; + else best_sign = 1; + } + } + + xi_sq = 0.0; + for (i = 0; i < 3; i++) { + tmp = displs[i] - best_sign * half_xi_chi_vec[which_crystal][best_i][i]; + xi_sq += tmp*tmp; + } + + if (xi_sq > 0.0) { + double xi = sqrt(xi_sq); + for (i = 0; i < 3; i++) + dxi[i] = (best_sign * half_xi_chi_vec[which_crystal][best_i][i] - + displs[i]) / xi; + } else dxi[0] = dxi[1] = dxi[2] = 0.0; +} + +/* ---------------------------------------------------------------------- + compare two neighbors I and J in sort data structure + called via qsort in post_force() method + is a static method so can't access sort data structure directly + return -1 if I < J, 0 if I = J, 1 if I > J + do comparison based on rsq distance +------------------------------------------------------------------------- */ + +int FixOrientBCC::compare(const void *pi, const void *pj) +{ + FixOrientBCC::Sort *ineigh = (FixOrientBCC::Sort *) pi; + FixOrientBCC::Sort *jneigh = (FixOrientBCC::Sort *) pj; + + if (ineigh->rsq < jneigh->rsq) return -1; + else if (ineigh->rsq > jneigh->rsq) return 1; + return 0; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double FixOrientBCC::memory_usage() +{ + double bytes = nmax * sizeof(Nbr); + bytes += 2*nmax * sizeof(double); + return bytes; +} diff --git a/src/MISC/fix_orient_bcc.h b/src/MISC/fix_orient_bcc.h new file mode 100644 index 0000000000..aecbaca793 --- /dev/null +++ b/src/MISC/fix_orient_bcc.h @@ -0,0 +1,115 @@ +/* -*- 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 FIX_CLASS + +FixStyle(orient/bcc,FixOrientBCC) + +#else + +#ifndef LMP_FIX_ORIENT_BCC_H +#define LMP_FIX_ORIENT_BCC_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixOrientBCC : public Fix { + public: + struct Nbr { // neighbor info for each owned and ghost atom + int n; // # of closest neighbors (up to 8) + tagint id[8]; // IDs of each neighbor + // if center atom is owned, these are local IDs + // if center atom is ghost, these are global IDs + double xismooth[8]; // distance weighting factor for each neighbors + double dxi[8][3]; // d order-parameter / dx for each neighbor + double duxi; // d Energy / d order-parameter for atom + }; + + struct Sort { // data structure for sorting to find 8 closest + int id; // local ID of neighbor atom + double rsq; // distance between center and neighbor atom + double delta[3]; // displacement between center and neighbor atom + double xismooth; // distance weighting factor + }; + + FixOrientBCC(class LAMMPS *, int, char **); + ~FixOrientBCC(); + int setmask(); + void init(); + void init_list(int, class NeighList *); + void setup(int); + void post_force(int); + void post_force_respa(int, int, int); + double compute_scalar(); + int pack_forward_comm(int, int *, double *, int, int *); + void unpack_forward_comm(int, int, double *); + double memory_usage(); + + private: + int me; + int nlevels_respa; + + int direction_of_motion; // 1 = center shrinks, 0 = center grows + int nstats; // stats output every this many steps + double a; // lattice parameter + double Vxi; // potential value + double uxif_low; // cut-off fraction, low order parameter + double uxif_high; // cut-off fraction, high order parameter + char *xifilename, *chifilename; // file names for 2 crystal orientations + + bool use_xismooth; + double Rxi[8][3],Rchi[8][3],half_xi_chi_vec[2][4][3]; + double xiid,xi0,xi1,xicutoffsq,cutsq,added_energy; + int half_bcc_nn; + + int nmax; // expose 2 per-atom quantities + double **order; // order param and normalized order param + + Nbr *nbr; + Sort *sort; + class NeighList *list; + + void find_best_ref(double *, int, double &, double *); + static int compare(const void *, const void *); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Fix orient/bcc file open failed + +The fix orient/bcc command could not open a specified file. + +E: Fix orient/bcc file read failed + +The fix orient/bcc command could not read the needed parameters from a +specified file. + +E: Fix orient/bcc found self twice + +The neighbor lists used by fix orient/bcc are messed up. If this +error occurs, it is likely a bug, so send an email to the +"developers"_http://lammps.sandia.gov/authors.html. + +*/