--- /dev/null 2004-02-18 08:26:44.000000000 -0700 +++ /home/sjplimp/lammps/src/fix_orient_fcc.cpp 2006-05-02 16:42:18.000000000 -0600 @@ -0,0 +1,552 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "stdlib.h" +#include "mpi.h" +#include "fix_orient_fcc.h" +#include "atom.h" +#include "update.h" +#include "respa.h" +#include "neighbor.h" +#include "comm.h" +#include "output.h" +#include "error.h" + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +/* ---------------------------------------------------------------------- */ + +FixOrientFCC::FixOrientFCC(int narg, char **arg) : Fix(narg, arg) +{ + MPI_Comm_rank(world,&me); + + if (narg != 11) error->all("Illegal fix orient/fcc command"); + + nstats = atoi(arg[3]); + direction_of_motion = atoi(arg[4]); + a = atof(arg[5]); + Vxi = atof(arg[6]); + uxif_low = atof(arg[7]); + uxif_high = atof(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("Illegal fix orient/fcc command"); + + // initializations + + PI = 4.0*atan(1.0); + half_fcc_nn = 6; + use_xismooth = false; + double xicutoff = 1.57; + xicutoffsq = xicutoff * xicutoff; + cutsq = 0.5 * a*a*xicutoffsq; + nmax = 0; + + // read xi and chi reference orientations from files + + if (me == 0) { + char line[512]; + char *result; + int count; + + FILE *infile = fopen(xifilename,"r"); + if (infile == NULL) error->one("Fix orient/fcc file open failed"); + for (int i = 0; i < 6; i++) { + result = fgets(line,512,infile); + if (!result) error->one("Read of fix orient/fcc file failed"); + count = sscanf(line,"%lg %lg %lg",&Rxi[i][0],&Rxi[i][1],&Rxi[i][2]); + if (count != 3) error->one("Read of fix orient/fcc file failed"); + } + fclose(infile); + + infile = fopen(chifilename,"r"); + if (infile == NULL) error->one("Fix orient/fcc file open failed"); + for (int i = 0; i < 6; i++) { + result = fgets(line,512,infile); + if (!result) error->one("Read of fix orient/fcc file failed"); + count = sscanf(line,"%lg %lg %lg",&Rchi[i][0],&Rchi[i][1],&Rchi[i][2]); + if (count != 3) error->one("Read of fix orient/fcc file 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 < 6; 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 12 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 < 6; 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 /= 12.0; + xi0 = uxif_low * xiid; + xi1 = uxif_high * xiid; +} + +/* ---------------------------------------------------------------------- */ + +FixOrientFCC::~FixOrientFCC() +{ + delete [] xifilename; + delete [] chifilename; + if (nmax) delete [] nbr; +} + +/* ---------------------------------------------------------------------- */ + +int FixOrientFCC::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= POST_FORCE_RESPA; + mask |= THERMO; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::init() +{ + if (use_xismooth) comm->maxforward_fix = MAX(comm->maxforward_fix,62); + else comm->maxforward_fix = MAX(comm->maxforward_fix,50); + + if (strcmp(update->integrate_style,"respa") == 0) + nlevels_respa = ((Respa *) update->integrate)->nlevels; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::setup() +{ + setup_flag = 1; + if (strcmp(update->integrate_style,"verlet") == 0) + post_force(1); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(1,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } + setup_flag = 0; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::post_force(int vflag) +{ + int i,j,k,m,n,nn,nsort,id_self; + int *neighs; + double edelta,added_energy,omega; + double dx,dy,dz,rsq,xismooth,xi_sq,duxi,duxi_other; + double dxi[3]; + double *dxiptr; + bool found_myself; + + int eflag = false; + if (thermo_print) { + if (setup_flag) eflag = true; + else if (output->next_thermo == update->ntimestep) eflag = true; + } + + // rebuild full neighbor list if normal neighbor list was just rebuilt + + if (neighbor->ago == 0) neighbor->build_full(); + + // set local ptrs + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + int *tag = atom->tag; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + + // insure nbr data structure is adequate size + + if (nall > nmax) { + if (nmax) delete [] nbr; + nbr = new Nbr[nall]; + nmax = nall; + } + + // loop over owned atoms and build Nbr data structure of neighbors + + if (eflag) added_energy = 0.0; + int count = 0; + int mincount = 1000000000; + int maxcount = 0; + + for (i = 0; i < nlocal; i++) { + neighs = neighbor->firstneigh_full[i]; + n = neighbor->numneigh_full[i]; + + if (n < mincount) mincount = n; + if (n > maxcount) { + if (maxcount) delete [] sort; + sort = new Sort[n]; + maxcount = n; + } + + // 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 (k = 0; k < n; k++) { + j = neighs[k]; + 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 - 2.0*rsq/(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 <= 12 + + if (nsort > 12) qsort(sort,nsort,sizeof(Sort),compare); + + // copy up to 12 nearest neighbors into nbr data structure + // operate on delta vector via find_best_ref() to compute dxi + + n = MIN(12,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; + + // compute potential derivative to xi + + if (xi_total < xi0) { + nbr[i].duxi = 0.0; + if (eflag) edelta = 0.0; + } else if (xi_total > xi1) { + nbr[i].duxi = 0.0; + if (eflag) edelta = Vxi; + } else { + omega = (0.5*PI)*(xi_total-xi0) / (xi1-xi0); + nbr[i].duxi = PI*Vxi*sin(2.0*omega) / (2.0*(xi1-xi0)); + if (eflag) edelta = Vxi*(1 - cos(2.0*omega)) / 2.0; + } + if (eflag) added_energy += edelta; + } + + if (maxcount) delete [] sort; + + // communicate to acquire nbr data for ghost atoms + + comm->comm_fix(this); + + // compute grain boundary force on each owned atom + // skip atoms not in group + + for (i = 0; i < nlocal; i++) { + 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("Fix orient/fcc 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]; + } + } + } + } + } + + // sum energy across procs + + if (eflag) + MPI_Allreduce(&added_energy,&total_added_e,1,MPI_DOUBLE,MPI_SUM,world); + + // 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 = 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 %d: %g atoms have %d neighbors\n", + update->ntimestep,atom->natoms,total); + fprintf(screen," neighs: min = %d, max = %d, ave = %g\n", + min,max,ave); + } + if (logfile) { + fprintf(logfile,"orient step %d: %g atoms have %d neighbors\n", + update->ntimestep,atom->natoms,total); + fprintf(logfile," neighs: min = %d, max = %d, ave = %g\n", + min,max,ave); + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +int FixOrientFCC::thermo_fields(int n, int *flags, char **keywords) +{ + if (n == 0) return 1; + flags[0] = 3; + strcpy(keywords[0],"Orient"); + return 1; +} + +/* ---------------------------------------------------------------------- */ + +int FixOrientFCC::thermo_compute(double *values) +{ + values[0] = total_added_e; + return 1; +} + +/* ---------------------------------------------------------------------- */ + +int FixOrientFCC::pack_comm(int n, int *list, double *buf, int *pbc_flags) +{ + int i,j,k,id,num; + int *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++] = static_cast (num); + buf[m++] = nbr[k].duxi; + + for (j = 0; j < num; j++) { + if (use_xismooth) buf[m++] = nbr[m].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++] = static_cast (id); + } + + m += (12-num) * 3; + if (use_xismooth) m += 12-num; + } + + if (use_xismooth) return 62; + return 50; +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::unpack_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++]); + } + + m += (12-num) * 3; + if (use_xismooth) m += 12-num; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixOrientFCC::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_fcc_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 FixOrientFCC::compare(const void *pi, const void *pj) +{ + FixOrientFCC::Sort *ineigh = (FixOrientFCC::Sort *) pi; + FixOrientFCC::Sort *jneigh = (FixOrientFCC::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 +------------------------------------------------------------------------- */ + +int FixOrientFCC::memory_usage() +{ + int bytes = atom->nmax * sizeof(Nbr); + return bytes; +} --- /dev/null 2004-02-18 08:26:44.000000000 -0700 +++ /home/sjplimp/lammps/src/fix_orient_fcc.h 2006-05-02 16:27:46.000000000 -0600 @@ -0,0 +1,77 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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. +------------------------------------------------------------------------- */ + +#ifndef FIX_ORIENT_FCC_H +#define FIX_ORIENT_FCC_H + +#include "fix.h" + +class FixOrientFCC : public Fix { + public: + struct Nbr { // neighbor info for each owned and ghost atom + int n; // # of closest neighbors (up to 12) + int id[12]; // IDs of each neighbor + // if center atom is owned, these are local IDs + // if center atom is ghost, these are global IDs + double xismooth[12]; // distance weighting factor for each neighbors + double dxi[12][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 12 closest + int id; // 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 + }; + + FixOrientFCC(int, char **); + ~FixOrientFCC(); + int setmask(); + void init(); + void setup(); + void post_force(int); + void post_force_respa(int, int, int); + int pack_comm(int, int *, double *, int *); + void unpack_comm(int, int, double *); + int thermo_fields(int, int *, char **); + int thermo_compute(double *); + int memory_usage(); + + private: + int me; + double PI; + 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; + int setup_flag; + double Rxi[12][3],Rchi[12][3],half_xi_chi_vec[2][6][3]; + double xiid,xi0,xi1,xicutoffsq,cutsq,total_added_e; + int half_fcc_nn,nmax; + + Nbr *nbr; + Sort *sort; + + void find_best_ref(double *, int, double &, double *); + static int compare(const void *, const void *); +}; + +#endif --- /home/sjplimp/oldlammps/src/style.h 2006-05-02 16:48:43.000000000 -0600 +++ /home/sjplimp/lammps/src/style.h 2006-05-02 16:28:39.000000000 -0600 @@ -80,11 +80,15 @@ #ifdef DumpInclude #include "dump_atom.h" #include "dump_custom.h" + //#include "dump_dcd.h" + //#include "dump_xtc.h" #endif #ifdef DumpClass DumpStyle(atom,DumpAtom) DumpStyle(custom,DumpCustom) + //DumpStyle(dcd,DumpDCD) + //DumpStyle(xtc,DumpXTC) #endif #ifdef FixInclude @@ -108,7 +112,7 @@ #include "fix_nvt.h" #include "fix_plane_force.h" #include "fix_print.h" - //#include "fix_orient_fcc.h" +#include "fix_orient_fcc.h" #include "fix_rdf.h" #include "fix_respa.h" #include "fix_rigid.h" @@ -147,7 +151,7 @@ FixStyle(npt,FixNPT) FixStyle(nve,FixNVE) FixStyle(nvt,FixNVT) - //FixStyle(orient/fcc,FixOrientFCC) +FixStyle(orient/fcc,FixOrientFCC) FixStyle(print,FixPrint) FixStyle(planeforce,FixPlaneForce) FixStyle(rdf,FixRDF) @@ -209,7 +213,6 @@ #include "pair_eam.h" #include "pair_eam_alloy.h" #include "pair_eam_fs.h" -#include "pair_hybrid.h" #include "pair_lj_cut.h" #include "pair_lj_cut_coul_cut.h" #include "pair_lj_cut_coul_debye.h" @@ -219,6 +222,8 @@ #include "pair_soft.h" #include "pair_table.h" #include "pair_yukawa.h" +#include "pair_hybrid.h" + //#include "pair_hybrid2.h" #endif #ifdef PairClass @@ -227,7 +232,6 @@ PairStyle(eam,PairEAM) PairStyle(eam/alloy,PairEAMAlloy) PairStyle(eam/fs,PairEAMFS) -PairStyle(hybrid,PairHybrid) PairStyle(lj/cut,PairLJCut) PairStyle(lj/cut/coul/cut,PairLJCutCoulCut) PairStyle(lj/cut/coul/debye,PairLJCutCoulDebye) @@ -237,6 +241,8 @@ PairStyle(soft,PairSoft) PairStyle(table,PairTable) PairStyle(yukawa,PairYukawa) +PairStyle(hybrid,PairHybrid) + //PairStyle(hybrid2,PairHybrid2) #endif #ifdef RegionInclude --- /dev/null 2004-02-18 08:26:44.000000000 -0700 +++ /home/sjplimp/lammps/doc/fix_orient_fcc.txt 2006-05-02 16:45:28.000000000 -0600 @@ -0,0 +1,135 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://www.cs.sandia.gov/~sjplimp/lammps.html) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +fix orient/fcc command :h3 + +fix ID group-ID orient/fcc nstats dir alat dE cutlo cuthi file0 file1 :pre + +ID, group-ID are documented in "fix"_fix.html command +nstats = print stats every this many steps, 0 = never +dir = 0/1 for which crystal is used as reference +alat = fcc cubic lattice constant (distance units) +dE = energy added to each atom (energy units) +cutlo,cuthi = values between 0.0 and 1.0, cutlo < cuthi +file0,file1 = files that specify orientation of each grain :ul + +[Examples:] + +fix gb all orient/fcc 0 1 4.032008 0.001 0.25 0.75 xi.vec chi.vec :pre + +[Description:] + +The fix applies an orientation-dependent force to atoms near a planar +grain boundary which can be used to induce grain boundary migration +(in the direction perpendicular to the grain boundary plane). The +motivation and explanation of this force and its application are +described in "(Janssens)"_#Janssens. The force is only applied to +atoms in the fix group. + +The basic idea is that atoms in one grain (on one side of the +boundary) have a potential energy dE added to them. Atoms in the +other grain have 0.0 potential energy added. Atoms near the boundary +(whose neighbor environment is intermediate between the two grain +orientations) have an energy between 0.0 and dE added. This creates +an effective driving force to reduce the potential energy of atoms +near the boundary by pushing them towards one of the grain +orientations. For dir = 1 and dE > 0, the boundary will thus move so +that the grain described by file0 grows and the grain described by +file1 shrinks. Thus this fix is designed for simulations of two-grain +systems, either with one grain boundary and free surfaces parallel to +the boundary, or a system with periodic boundary conditions and two +equal and opposite grain boundaries. In either case, the entire +system can displace during the simulation, and such motion should be +accounted for in measuring the grain boundary velocity. + +The potential energy added to atom I is given by these formulas + +:c,image(Eqs/fix_orient_fcc.jpg) + +which are fully explained in "(Janssens)"_#Janssens. The order +parameter Xi for atom I in equation (1) is a sum over the 12 nearest +neighbors of atom I. Rj is the vector from atom I to its neighbor J, +and RIj is a vector in the reference (perfect) crystal. That is, if +dir = 0/1, then RIj is a vector to an atom coord from file 0/1. +Equation (2) gives the expected value of the order parameter XiIJ in +the other grain. Hi and lo cutoffs are defined in equations (3) and +(4), using the input parameters {cutlo} and {cuthi} as threshholds to +avoid adding grain boundary energy when the deviation in the order +parameter from 0 or 1 is small (e.g. due to thermal fluctuations in a +perfect crystal). The added potential energy Ui for atom I is given +in equation (6) where it is interpolated between 0 and dE using the +two threshhold Xi values and the Wi value of equation (5). + +The derivative of this energy expression gives the force on each atom +which thus depends on the orientation of its neighbors relative to the +2 grain orientations. Only atoms near the grain boundary feel a net +force which tends to drive them to one of the two grain orientations. + +In equation (1), the reference vector used for each neigbbor is the +reference vector closest to the actual neighbor position. This means +it is possible two different neighbors will use the same reference +vector. In such cases, the atom in question is far from a perfect +orientation and will likely receive the full dE addition, so the +effect of duplicate reference vector usage is small. + +The {dir} parameter determines which grain wants to grow at the +expense of the other. A value of 0 means the first grain will shrink; +a value of 1 means it will grow. This assumes that {dE} is positive. +The reverse will be true if {dE} is negative. + +The {alat} parameter is the cubic lattice constant for the fcc +material and is only used to compute a cutoff distance of 1.57 * alat +/ sqrt(2) for finding the 12 nearest neighbors of each atom (which +should be valid for an fcc crystal). A longer/shorter cutoff can be +imposed by adjusting {alat}. If a particular atom has less than 12 +neighbors within the cutoff, the order parameter of equation (1) is +effectively multiplied by 12 divided by the actual number of neighbors +within the cutoff. + +The {dE} parameter is the maximum amount of additional energy added to +each atom in the grain which wants to shrink. + +The {cutlo} and {cuthi} parameters are used to reduce the force added +to bulk atoms in each grain far away from the boundary. An atom in +the bulk surrounded by neighbors at the ideal grain orientation would +compute an order parameter of 0 or 1 and have no force added. +However, thermal vibrations in the solid will cause the order +parameters to be greater than 0 or less than 1. The cutoff parameters +mask this effect, allowing forces to only be added to atoms with +order-parameters between the cutoff values. + +{File0} and {file1} are filenames for the two grains which each +contain 6 vectors (6 lines with 3 values per line) which specify the +grain orientations. Each vector is a displacement from a central atom +(0,0,0) to a nearest neighbor atom in an fcc lattice at the proper +orientation. The vector lengths should all be identical since an fcc +lattice has a coordination number of 12. Only 6 are listed due to +symmetry, so the list must include one from each pair of +equal-and-opposite neighbors. + +This fix supports the "fix_modify"_fix_modify.html options for +{thermo} and {energy}. The former will print the contribution the fix +makes to the energy of the system when thermodynamics is printed. The +latter will add this contribution to the total potential energy +(PotEng). + +[Restrictions:] + +This fix should only be used with fcc lattices. + +[Related commands:] + +"fix_modify"_fix_modify.html + +[Default:] none + +:line + +:link(Janssens) +[(Janssens)] Janssens, Olmsted, Holm, Foiles, Plimpton, Derlet, Nature +Materials, 5, 124-127 (2006). --- /dev/null 2004-02-18 08:26:44.000000000 -0700 +++ /home/sjplimp/lammps/doc/fix_orient_fcc.html 2006-05-02 16:45:30.000000000 -0600 @@ -0,0 +1,141 @@ + +
LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
+ + + + + + +
+ +

fix orient/fcc command +

+
fix ID group-ID orient/fcc nstats dir alat dE cutlo cuthi file0 file1 
+
+
  • ID, group-ID are documented in fix command +
  • nstats = print stats every this many steps, 0 = never +
  • dir = 0/1 for which crystal is used as reference +
  • alat = fcc cubic lattice constant (distance units) +
  • dE = energy added to each atom (energy units) +
  • cutlo,cuthi = values between 0.0 and 1.0, cutlo < cuthi +
  • file0,file1 = files that specify orientation of each grain +
+

Examples: +

+
fix gb all orient/fcc 0 1 4.032008 0.001 0.25 0.75 xi.vec chi.vec 
+
+

Description: +

+

The fix applies an orientation-dependent force to atoms near a planar +grain boundary which can be used to induce grain boundary migration +(in the direction perpendicular to the grain boundary plane). The +motivation and explanation of this force and its application are +described in (Janssens). The force is only applied to +atoms in the fix group. +

+

The basic idea is that atoms in one grain (on one side of the +boundary) have a potential energy dE added to them. Atoms in the +other grain have 0.0 potential energy added. Atoms near the boundary +(whose neighbor environment is intermediate between the two grain +orientations) have an energy between 0.0 and dE added. This creates +an effective driving force to reduce the potential energy of atoms +near the boundary by pushing them towards one of the grain +orientations. For dir = 1 and dE > 0, the boundary will thus move so +that the grain described by file0 grows and the grain described by +file1 shrinks. Thus this fix is designed for simulations of two-grain +systems, either with one grain boundary and free surfaces parallel to +the boundary, or a system with periodic boundary conditions and two +equal and opposite grain boundaries. In either case, the entire +system can displace during the simulation, and such motion should be +accounted for in measuring the grain boundary velocity. +

+

The potential energy added to atom I is given by these formulas +

+
+
+

which are fully explained in (Janssens). The order +parameter Xi for atom I in equation (1) is a sum over the 12 nearest +neighbors of atom I. Rj is the vector from atom I to its neighbor J, +and RIj is a vector in the reference (perfect) crystal. That is, if +dir = 0/1, then RIj is a vector to an atom coord from file 0/1. +Equation (2) gives the expected value of the order parameter XiIJ in +the other grain. Hi and lo cutoffs are defined in equations (3) and +(4), using the input parameters cutlo and cuthi as threshholds to +avoid adding grain boundary energy when the deviation in the order +parameter from 0 or 1 is small (e.g. due to thermal fluctuations in a +perfect crystal). The added potential energy Ui for atom I is given +in equation (6) where it is interpolated between 0 and dE using the +two threshhold Xi values and the Wi value of equation (5). +

+

The derivative of this energy expression gives the force on each atom +which thus depends on the orientation of its neighbors relative to the +2 grain orientations. Only atoms near the grain boundary feel a net +force which tends to drive them to one of the two grain orientations. +

+

In equation (1), the reference vector used for each neigbbor is the +reference vector closest to the actual neighbor position. This means +it is possible two different neighbors will use the same reference +vector. In such cases, the atom in question is far from a perfect +orientation and will likely receive the full dE addition, so the +effect of duplicate reference vector usage is small. +

+

The dir parameter determines which grain wants to grow at the +expense of the other. A value of 0 means the first grain will shrink; +a value of 1 means it will grow. This assumes that dE is positive. +The reverse will be true if dE is negative. +

+

The alat parameter is the cubic lattice constant for the fcc +material and is only used to compute a cutoff distance of 1.57 * alat +/ sqrt(2) for finding the 12 nearest neighbors of each atom (which +should be valid for an fcc crystal). A longer/shorter cutoff can be +imposed by adjusting alat. If a particular atom has less than 12 +neighbors within the cutoff, the order parameter of equation (1) is +effectively multiplied by 12 divided by the actual number of neighbors +within the cutoff. +

+

The dE parameter is the maximum amount of additional energy added to +each atom in the grain which wants to shrink. +

+

The cutlo and cuthi parameters are used to reduce the force added +to bulk atoms in each grain far away from the boundary. An atom in +the bulk surrounded by neighbors at the ideal grain orientation would +compute an order parameter of 0 or 1 and have no force added. +However, thermal vibrations in the solid will cause the order +parameters to be greater than 0 or less than 1. The cutoff parameters +mask this effect, allowing forces to only be added to atoms with +order-parameters between the cutoff values. +

+

File0 and file1 are filenames for the two grains which each +contain 6 vectors (6 lines with 3 values per line) which specify the +grain orientations. Each vector is a displacement from a central atom +(0,0,0) to a nearest neighbor atom in an fcc lattice at the proper +orientation. The vector lengths should all be identical since an fcc +lattice has a coordination number of 12. Only 6 are listed due to +symmetry, so the list must include one from each pair of +equal-and-opposite neighbors. +

+

This fix supports the fix_modify options for +thermo and energy. The former will print the contribution the fix +makes to the energy of the system when thermodynamics is printed. The +latter will add this contribution to the total potential energy +(PotEng). +

+

Restrictions: +

+

This fix should only be used with fcc lattices. +

+

Related commands: +

+

fix_modify +

+

Default: none +

+
+ + + +

(Janssens) Janssens, Olmsted, Holm, Foiles, Plimpton, Derlet, Nature +Materials, 5, 124-127 (2006). +

+ --- /dev/null 2004-02-18 08:26:44.000000000 -0700 +++ /home/sjplimp/lammps/doc/Eqs/fix_orient_fcc.jpg 2006-05-02 16:37:03.000000000 -0600 @@ -0,0 +1,27 @@ +JFIFHHC  !"$"$ Q +!1"AQa#2Vq$BՁ%3SUrRbcu&57Cegs?gƸme"|foLI[vAx:' XK5ڵm=@`f;ek\|΁ysp~u_/ūeyߎ,qKO <1psɲwh;YR5JbwC s Ҏ8~WǪr sm0M鹝8Av%V񵖬V!hf|N?Ж?B""cxIQ?2 V&E+Nɯ-XDM%P.;WVR+tCfXǴGh?)W-* X菌=XGQ}kp8Qڡv^fc])ըZDlJfgw!`p\)S壹1!Q,i6kIx+Y9>j<}mֹW8s\?cʹEO97q1}̶ FG1GICcxv0՝v<:^s\H ZnW#s;i2%/|~OigLݸKNnck}Gd.nֶwFgnZ3i5`QXSxVN\6{RcnFm7}:n hQ_.~ž&\8$KAe7옩*Jc\>7݀Ӽ$%g:2Kohhp-~`/"(2AR(e5:iCC? k[YnEn"q3 xNj$h=x$ydg>2#"(Sb(Mh n)ȱ\.cw ✰{ؖly^egrNEkcًYv|6:F??-~^n>T!kmW0˲_\Z=ɳ7yشRrTOob6??bxcy>W=p$˙*N`pbZǟ?WXJVۭYY sZz6O]l[CZ#zib\vH:WUv.%s ף:#o df˗>c ^,tykW5l_Y whۢ󭯟*ʍy|a?68~lMx.ޜlgw;@O깢,S'¬nլ\qbyts瑢w{];N_ms4I48*b*1aUK~5.!B7={M -da{J6"وr#"&gf83}C,C芷ٜƻ5u]h8CNۣ>cZq1 XjXh}5DDDDDDDDDDD\ QZ K{K\uɠ54""""""",5{D3vE(ՙM>`HdoK% >b;;.*orJkl@]My|:bԤvq=~'8ZDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDPb%2tot>5VV;TCfl'x+y +sŹW7n>XORb nۚu<úRUᱎ$,Oc_RYq`.s˺7:Wo]#,K2~1]Yn'KZfJvǷ4ڈy 5!Ԟ *ԁMmf'@xB/Z)9"Fވ]˿ ].K29qRj䧫Z,} ŏ9Q9i]<5܊|Otooi#qXIrә.C=;Xܞ^f9>z +^FHE; s[$ cAfx*fX6$pkZ?2O lY͋6sנ3- ^ 5ZXWM;bpk{ +gÉ#Ztok3YVG!i=AFGX 2;Xp޺*Vy$gҚdPB,uw^1-{Ya20屑^`qsqc湧^6_eӚd=dkcn}Mȷ\0 Ɇ,lW^/,N]n"{ƑR4t=Q9'V>C__oJw gֳ/^w+oO,PBEKlO ~)X䖣O}H[c~@@ ·x䍆8Wԕ}t۶5O jy9CZ!)Oצw|g?|OV)!Agumk o֢""",28TʿVQ숇]IQsOwe8'vRYfꎦ2\Aszuh;ڛU}l|oI߷|8sI$}[rKJ=Ԩ[d;-2#3ۆ:oawk s,~Zn:A #ss'AZ9XOO+GA'cGi,㌧ϝh +Tݽ[{ekC, ikwZ$z (Lt8dT`s\I#^9q׍SQRf= ZbY{l Lqr}__9˼= ⎦֟2nCl$Jƻ{ + P%(. 'M=wg4pkI-`>{ e>vTDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDEkױp\-\Uը t$ 285?]1{Zѷ8?5p2&ipo@B+NG.q?`=&4#+V\UḼ5!pl8_ \Hv5ώ\b48ثY +r wN'g6G߾' NY=/drn(Z}qhfVnGC7J7:yXԆNvy"0b)'{uU5lQ3gA{>u[6kqLr8G2$]X\Om!sz뽍Oc%L}FYad,s%i`h-jǰpՖZcLF3d𲗶NΌhOϠɳufYg-'1ԇuiۛص<>%OWS2o6|zr69PlOs44\;nn ҵ)C>*&綼apy/>@>;ȕ7eQWcs oW=8u ׿VUy۳Mo6. P RA@sυMĪԡ5]I~CnHuGp\@O4 Ë]eC~#-=(w鑍+X?FDDDDYq_Z$q\s%iMab77$]^5Gs7̹%FbYaoаbsRFDzۻ3嫷=ryjJNxnR{VEX`FCO%_|82qlX61;>s&3;B* c^oKoOz^zZCGrTsAB&=a +WrL&?bE4oc)cx|r1g5kxQz[{ raesָ71h;yjJ<%'6O>S:191&@ZXC\}`}cZJ1 =:c=:NΆdU%q6ψ%d 3}9}ču"G-&F:'YVZִonvn֗ +Kٱj! ̕ٯXg^ ZU""""KwDDDD]SYj:wl}.$dL#Ƃt< 5km`}DDDDDDDDDDDECz} dQ:Y4\C#`""""/-eP2E5 Ø_(@}{0mdJZ84qu{7f{'*/*T|g)ogTAK@$l$U٬9_So| vegk~z3Y+7#n!Ν<,gC's +TϑUv_=(ᠽd :ҭ5Ǹg9wTqk[b..pw9K|'1yF>,0ղUgw9װ45I8yjd+7edMQtShI'zC~Әrjw8_01K&;z8.P~X":}ѝsHӕYsCW32,ׯ3+L7 &-q ^CaffE G!asJ3Iw+x񿗲Pe! X4CvųK9O#J#dyYiX&Fݟ#>mOn̔`,f݁ 'W-w x?F4~l`U#Rjb7F^k|\ h)oمդϋcf,F]ݍneHA3 qೞdD-HL-wQs.[#y ,,Md)е5^l0FB,u&mq&CVf" W;-٥\o}U٬7~3dޗdޞ7(8܎/f#悄L{ÃvڮL ;΋Վhu|Rck8l-n.罭poۨcZwz8XyK7>O%rl: }tb7s cK]"Mxߕ 032VCZ2w:X÷6@4/Pe%7zAcwt?}wLlu/a򏥖ƹW?zEvizk큣5rNUW=/1&6>\F? ˷yIhttbeb6y_^Noo N82~۶r9Fz:C_sɑ&M3m;`kk +Kٱj! ̕ٯXg^ ZU"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""/08A ~Z2#~ᤀHU|esݍۜvO'R'U>N${K9vx~s^5Oc2VkYtpa{ӧ?gz>eVid֭L,p%Tr3+@-/u%z-cԀJDDDDDDDDDDDDDEU&j5i,=6d0ϾqSOF{509_3LHk |VL9ldWk\G#X;i׍SQyvqF|[^$mW/m_칢0y S^TI2V'n5{]_!>\ֻ{>Z~:]>s[!}=?={޼c?'Au^7<.k'~ėޝygq)Չdc74 Z Wlc?:Pkۏ 3+Fϧm|sPo7 Գ7=>FɮN.׶W?~)9pK? o?.gA? ?,. Odk ߊ\3R~OˆjYPO~)9pK? M+.wi];zӨ?:M;yބDDDDDDDDDDDDEzkqNk-6/2ա{,'Q94O'f#fͬ+&,Ecc`ӜO^΁+`~_Z^ֿ]j" 9=^jqt]Տ{õ܎޵ +b~!/[3uzb4f5ѿLsAϛ#6]/c0:4@c,zm2x; s?UܶVS [=Ib25Eqdb@IZV$tP7=瑱:8XN [0cIYlqkA>6B͙A++F9t>i%wӹNv֝E xroOOb(1bwcgM0inBc7B1natrtx<˺HtR4cP""""""""""""""(1zGֆ XQםevWFS\Ϻc"GImslҚ ˯] :+pr\-Vbߟ%O_ra%1>@ FK)!ј{4]3.R{1=I$O +O˼ДU=>jׂ(C^{8sDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDESx\Ӻ$V=Mda΍#IBDDDDDDDTy= +܅6׿jVea)x[ o E9aj2Ju|=t I)ughJ;Dml?Ejx*֖թXdY4 |JQtXu$;@G=y/٬.KK^rÇE-jύi/&6{F' b4ݰ]Hd{leĜ~CV4ok:xmpE7$DDDDDDDE +#vwsJ̯s׹GNuVկ1 1/AZ4`d{Z +XlU+Rxl?/c =99gy;>TDDDDDDE + F&FJ ](Z`e,3F66|u`U|.6:k$cv4@>Rq((GN*6(+cGkZ,1fR6/1z3!Vy"4@xi3$r7q0e2OaLɥc FAC;'ʰPc˔q\&1iX$lC}GGPjrڷ3l-LnF{ƶ07%;m{ + -!YGFc5>ZAbn +QvU{C&*my{֜=~KOu>R{*8ϥ3%Ͽb5s&G#j-G\Go̱}=[  }':PUWVH{)>2Da2c9-!q;F3&Ve298\c>fi$dr7`<8;%ا_c217ԯ64&o]rr4K3e;~ڠi,zp.h=@^O]h/JԉԭA^~bgtG# 2W~گn7Ƥ;l%-kf21ߪs۳s rNyHx|Q9?~AeESy, U͈YjM1ln[6HƤR,ddl7Fuڢ48ٲ +n{6N@H 6\(җՂMR7A|""f_8+ݹph:oI .Xr4J_V 7Hމ +U-SG,85h.$$錱vIivןЏo +z"""(+ɒF7v}hX I4O rL>Jc[m,PEFv ieީ'xfofh| HDDUܒՊx[jν_|us>=πV3گAб;%з#ܯG |q͍ lK=_$^/ʛlllu&M$yur1Bf*5:7hK`f7Mv:YfXGWHב.^ZBs|SMv/ɼ3z~%O7uOOľ)&??g=??V>&t.p> D;\>/ɼ3z~%O7uOOľ)&?뿋An +y>-1tS٘9ݝFq<%|b,/;Y>$'7펌 ׊i>}|H9ys2vv<{#⫢O1~kHSrI{Mb75bU]VE Hy 7awl惧 ck]sKtNL`@s)qˑ&)ޖOXvYY˝ςO0017ci}ZpOBR[{YcѼ pk؏YYS.ץ}Cϫ(s>Ws\ܥjK־K8-!Q!T</ wF5ntc!|-kix.ié<3UJMt [$q|k[*8pɈ|y:؛E|?a*/I DUys tZ\To 5N*-c'vCk} ͍#,roy/,\ܼ?- ~^Z-;[Zw0[,ӆ33{i` UV'uFuͳӏ0EwSix/3ck82*9io~m^| yb2|FĢ 8eƈ rcK{xo^|9*86./.׶o4=(CČ|cDB]Ve+6黵i偯hmD~hF%n.^^l{'/ bqؗJKRz^c2@vqh%U͸Kw!Uխ׵-i޶,D8KNATvo\Ye6?\2N%_U;N:r45*6v?d翙3ؑ"8̅ؖΆ $XqxbnFħ xog;CއT,~+fz5#ND:[>6N`[$7<;et i]?nFV_>/|/)qX2eR1~&6GWGWMOSl;七Fq!cİHdHd$e$vi1{# Qc:Y dޭi\$/6~$l&CRS;ECFi=1(2t;~]e0Zu>?t^C:mvq؊I]^Dz g[׶R^;9E$d~msXb=vA І醵GV`gMXbT&bġ}Iy|;8yұDDDDDDDDv,O>o;Y9rf{3=kq27`c=q+lq ~//ZF^C!YcVѷ0$ Wb+*od1v. y,=֍k2M̼u{$l^C^׳@:DDDDDDX>9&,LU\>+$ޫ}(${Ah[c7[||8t->w'^=~>rF+^t[} ;Wr)Y0sC_#-$n藑kT"^6G9Yဴ!t$kdusE63R֒d׆844Ҳ\$`.󡳣)!mfMC nKt=짢""""""""*I_ cz9 +5Ӄ7hwÿzp=0uhczq>jC΅FFG!>wknFSVuΰ!JGo|nӾ'scMbs=;7O #z>X|^2)$7c5;Eo@+$DDDDDD^OCrx.YeZߔ6!ݘg34;1kOW;[:x""ܞВ#ֈLmx`=#<>NרҊ(/dZ y.O벰\Z^Ge.Z'yrf/^ \@wfVsYo7~5He"!bƎojز?pѐ{%ˡA#_Ӹ$M~~H\Ǐ~ӆ'޽zsDZv7NbXFVo^YR oF:$k%!BKkY|,d֬/K t-k\vt@+ː\oVjFb_ՀK7[lwQ뭺go=Fu*vknYiq6÷'[ <.fUƭx*+أi$5NI$~䒻Q \ No newline at end of file