half-bond lists and no bond migration for local hyper

using namespace LAMMPS_NS;

using namespace FixConst;

#define DELTA 16384

#define DELTABOND 16384

#define VECLEN 5

// NOTE: count/output # of timesteps on which bias is non-zero

// NOTE: should there be a virial contribution from boosted bond?

// NOTE: allow newton off?  see Note in pre_reverse()

// possible enhancements

//   should there be a virial contribution from boosted bond?

//   allow newton off?  see Note in pre_reverse()

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

{

  int mask = 0;

  mask |= PRE_NEIGHBOR;

  mask |= PRE_FORCE;

  mask |= PRE_REVERSE;

  mask |= THERMO_ENERGY;

  return mask;

  if (force->newton_pair == 0)

    error->all(FLERR,"Hyper global requires newton pair on");

  if (atom->molecular && me == 0)

    error->warning(FLERR,"Hyper global for molecular systems "

                   "requires care in defining hyperdynamic bonds");

  dt = update->dt;

  // need an occasional half neighbor list

void FixHyperGlobal::pre_neighbor()

{

  int m,iold,jold,ilocal,jlocal;

  int i,m,iold,jold,ilocal,jlocal;

  double distsq;

  // reset local IDs for owned bond atoms, since atoms have migrated

  // uses xold and tagold from when bonds were created

  // reset local indices for owned bond atoms, since atoms have migrated

  // must be done after ghost atoms are setup via comm->borders()

  // first time this is done for a particular I or J atom:

  //   use tagold and xold from when bonds were created

  //   atom->map() finds atom ID if it exists, owned index if possible

  //   closest current I or J atoms to old I may now be ghost atoms

  //   closest_image() returns the ghost atom index in that case

  // also compute max drift of any atom in a bond

  //   drift = displacement from quenched coord while event has not yet occured

  for (i = 0; i < nall_old; i++) old2now[i] = -1;

  double **x = atom->x;

  for (m = 0; m < nblocal; m++) {

    iold = blist[m].iold;

    jold = blist[m].jold;

    ilocal = old2now[iold];

    jlocal = old2now[jold];

    if (ilocal < 0) {

      ilocal = atom->map(tagold[iold]);

      ilocal = domain->closest_image(xold[iold],ilocal);

      if (ilocal < 0)

        error->one(FLERR,"Fix hyper/global bond atom not found");

      old2now[iold] = ilocal;

      distsq = MathExtra::distsq3(x[ilocal],xold[iold]);

      maxdriftsq = MAX(distsq,maxdriftsq);

    }

    if (jlocal < 0) {

      jlocal = atom->map(tagold[jold]);

      jlocal = domain->closest_image(xold[iold],jlocal);   // closest to iold

      if (jlocal < 0)

        error->one(FLERR,"Fix hyper/global bond atom not found");

      old2now[jold] = jlocal;

      distsq = MathExtra::distsq3(x[jlocal],xold[jold]);

      maxdriftsq = MAX(distsq,maxdriftsq);

    }

    blist[m].i = ilocal;

    blist[m].j = jlocal;

  }

  /* old way - nblocal loop is re-doing index-find calculation

  // NOTE: drift may not include J atoms moving (if not themselves bond owners)

  int flag = 0;

  for (m = 0; m < nblocal; m++) {

  }

  if (flag) error->one(FLERR,"Fix hyper/global bond atom not found");

  */

}

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

{

  int i,j,m,imax,jmax;

  double delx,dely,delz;

  double r,r0,estrain,rmax,r0max,emax,dt_boost;

  double vbias,fbias,fbiasr;

  double r,r0,estrain,rmax,r0max,dt_boost;

  double ebias,vbias,fbias,fbiasr;

  // compute current strain of each owned bond

  // emax = maximum strain of any bond I own

  // eabs_max = maximum absolute value of strain of any bond I own

  // imax,jmax = local indices of my 2 atoms in that bond

  // rmax,r0max = current and relaxed lengths of that bond

  double **x = atom->x;

  emax = 0.0;

  double estrain_maxabs = 0.0;

  for (m = 0; m < nblocal; m++) {

    i = blist[m].i;

    r0 = blist[m].r0;

    estrain = fabs(r-r0) / r0;

    if (estrain > emax) {

      emax = estrain;

    if (estrain > estrain_maxabs) {

      estrain_maxabs = estrain;

      rmax = r;

      r0max = r0;

      imax = i;

  // finds max strain and what proc owns it

  // owner = proc that owns that bond

  pairme.value = emax;

  pairme.value = estrain_maxabs;

  pairme.proc = me;

  MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MAXLOC,world);

  owner = pairall.proc;

    return;

  }

  // I own the bond with max strain

  // compute Vbias and apply force to atoms imax,jmax

  // NOTE: logic would need to be different for newton off

  // I own the bond with max absolute value of strain

  // compute bias force on atoms imax,jmax if strain < q, else zero

  // Ebias = current strain = (r-r0) / r0

  // Vbias = bias potential = Vmax (1 - Ebias^2/q^2)

  // Fbias = bias force as function of strain

  //       = -dVbias/dEbias = 2 Vmax Ebias / q^2

  // Fix = x component of force on atom I

  //     = Fbias dEbias/dr dr/dxi, dEbias/dr = 1/r0, dr/dxi = delx/r

  // dt_boost = time boost factor = exp(Vbias/kT)

  // NOTE: logic here would need to be different for newton off

  double **f = atom->f;

  vbias = fbias = 0.0;

  dt_boost = 1.0;

  if (emax < qfactor) {

    vbias = vmax * (1.0 - emax*emax*invqfactorsq);

    fbias = 2.0 * vmax * emax / (qfactor*qfactor * r0max);

  if (estrain_maxabs < qfactor) {

    //ebias = (rmax-r0max) / r0max;

    ebias = fabs(rmax-r0max) / r0max;

    vbias = vmax * (1.0 - ebias*ebias*invqfactorsq);

    fbias = 2.0 * vmax * ebias * invqfactorsq;

    dt_boost = exp(beta*vbias);

    delx = x[imax][0] - x[jmax][0];

    dely = x[imax][1] - x[jmax][1];

    delz = x[imax][2] - x[jmax][2];

    fbiasr = fbias / rmax;

    fbiasr = fbias / r0max / rmax;

    f[imax][0] += delx*fbiasr;

    f[imax][1] += dely*fbiasr;

    f[imax][2] += delz*fbiasr;

    f[jmax][0] -= delx*fbiasr;

    f[jmax][1] -= dely*fbiasr;

    f[jmax][2] -= delz*fbiasr;

  } else nobias++;

  // output quantities

  outvec[0] = vbias;

  outvec[1] = dt_boost;

  outvec[2] = emax;

  outvec[2] = ebias;

  outvec[3] = atom->tag[imax];

  outvec[4] = atom->tag[jmax];

  if (atom->nmax > maxold) {

    memory->destroy(xold);

    memory->destroy(tagold);

    memory->destroy(old2now);

    maxold = atom->nmax;

    memory->create(xold,maxold,3,"hyper/global:xold");

    memory->create(tagold,maxold,"hyper/global:tagold");

    memory->create(old2now,maxold,"hyper/global:old2now");

  }

  tagint *tag = atom->tag;

  int nall = atom->nlocal + atom->nghost;

  nall_old = nall;

  for (i = 0; i < nall; i++) {

    xold[i][0] = x[i][0];

void FixHyperGlobal::grow_bond()

{

  // NOTE: could add int arg to do initial large alloc:

  // maxbond = maxbond/DELTA * DELTA; maxbond += DELTA;

  maxbond += DELTA;

  if (maxbond < 0 || maxbond > MAXSMALLINT)

    error->one(FLERR,"Fix hyper/local per-processor bond count is too big");

  if (maxbond + DELTABOND > MAXSMALLINT)

    error->one(FLERR,"Fix hyper/global bond count is too big");

  maxbond += DELTABOND;

  blist = (OneBond *)

    memory->srealloc(blist,maxbond*sizeof(OneBond),"hyper/local:blist");

    memory->srealloc(blist,maxbond*sizeof(OneBond),"hyper/global:blist");

}

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

  // 11 vector outputs returned for i = 0-10

  // i = 0 = boost factor on this step

  // i = 1 = max strain of any bond on this step

  // i = 1 = max strain of any bond on this step (positive or negative)

  // i = 2 = ID of atom I in max-strain bond on this step

  // i = 3 = ID of atom J in max-strain bond on this step

  // i = 4 = ave bonds/atom on this step

  if (i == 3) return outvec[4];

  if (i == 4) {

    int allbonds;     // NOTE: bigint?

    MPI_Allreduce(&nblocal,&allbonds,1,MPI_INT,MPI_SUM,world);

    bigint mybonds = nblocal;

    bigint allbonds;

    MPI_Allreduce(&mybonds,&allbonds,1,MPI_LMP_BIGINT,MPI_SUM,world);

    return 2.0*allbonds/atom->natoms;

  }

  // coords and IDs of owned+ghost atoms when bonds were formed

  // persists on a proc from one event until the next

  int nall_old;                 // nlocal+nghost for old atoms

  int maxold;                   // allocated size of old atoms

  double **xold;                // coords of atoms when bonds were formed

  tagint *tagold;               // IDs of atoms when bonds were formed

  tagint *tagold;               // IDs of atoms when bonds were forme

  int *old2now;                 // o2n[i] = current local index of old atom I

  // MPI data struct for finding bond with max strain via Allreduce

  int pack_reverse_comm(int, int, double *);

  void unpack_reverse_comm(int, int *, double *);

  void grow_arrays(int);

  void copy_arrays(int, int, int);

  int pack_exchange(int, double *);

  int unpack_exchange(int, double *);

  double memory_usage();

  // extra methods visible to callers

  double cutbond,qfactor,vmax,tequil,dcut;

  double alpha_user;         // timescale to apply boostostat (time units)

  double alpha;              // unitless dt/alpha_user

  double boosttarget;        // target value of boost

  int histoflag;

  int lostbond,lostbond_partner;

  double lostbond_coeff;

  double boost_target;       // target value of boost

  int checkbias,checkbias_every,checkbias_flag,checkbias_count;

  int checkcoeff,checkcoeff_every,checkcoeff_flag,checkcoeff_count;

  int setupflag;             // 1 during setup, 0 during run

  int firstflag;             // set for first time bond_build takes place

  int nostrainyet;           // 1 until maxstrain is first computed

  int nboost_running,nobias_running;

  int nbias_running,nobias_running;

  int nbondbuild;

  double time_bondbuild;

  bigint starttime;

  double sumboostcoeff;  // sum of aveboost at every timestep

  int allbonds;          // sum of bond count on this step

  double allboost;       // sum of boostcoeff on all bonds on this step

  double sumbiascoeff;   // sum of aveboost at every timestep

  bigint allbonds;       // sum of bond count on this step

  double allbias;        // sum of biascoeff on all bonds on this step

  int nnewbond;              // running tally of number of new bonds created

  int maxbondperatom;        // max # of bonds any atom ever has

  double mybias;

  double maxbondlen;         // cummulative max length of any bond

  double maxdriftsq;         // max distance any atom drifts from original pos

  double maxboostcoeff;      // cummulative max boost coeff for any bond

  double minboostcoeff;      // cummulative min boost coeff for any bond

  double maxbiascoeff;       // cummulative max bias coeff for any bond

  double minbiascoeff;       // cummulative min bias coeff for any bond

  double rmaxever,rmaxeverbig;

  int ghost_toofar;

  // extra timers

  class NeighList *listfull;   // full neigh list up to Dcut distance

  class NeighList *listhalf;   // half neigh list up to pair distance

                               // both created only when bonds are rebuilt

  //double timefirst,timesecond,timethird,timefourth;

  //double timefifth,timesixth,timeseventh,timetotal;

  // list of my owned bonds

  // persists on a proc from one event until the next

  // data structs for per-atom and per-bond info

  // all of these are for current owned and ghost atoms

  // except list and old2now are for atom indices at time of last bond build

  struct OneBond {             // single IJ bond, atom I is owner

    int i,j;                   // current local indices of 2 bond atoms

    int iold,jold;             // local indices when bonds were formed

    double r0;                 // relaxed bond length

    double biascoeff;          // biasing coefficient = prefactor Cij

  };

  class NeighList *list;       // full neigh list up to Dcut distance

                               // created only when bonds are rebuilt

  struct OneBond *blist;       // list of owned bonds

  int nblocal;                 // # of owned bonds

  int maxbond;                 // allocated size of blist

  int *old2now;                // o2n[i] = current local index of old atom i

                               //   stored for old owned and ghost atoms

                               //   I = old index when bonds were last created

                               //   old indices are stored in old neighbor list

  // old data from last timestep bonds were formed

  // persists on a proc from one event until the next

  // first set of vectors are maxlocal in length

  // second set of vectors are maxall in length

  double **xold;               // coords of owned+ghost atoms when bonds created

  tagint *tagold;              // global IDs of owned+ghost atoms when b created

  int nlocal_old;               // nlocal for old atoms

  int nall_old;                 // nlocal+nghost for old atoms

  int maxlocal;                 // allocated size of old local atom vecs

  int maxall;                   // allocated size of old all atom vecs

  int maxold;                  // allocated size of old2now

  int maxbond;                 // allocated size of bonds

  int old_nall;                // nlocal+nghost when old2now was last setup

  int *numbond;                 // # of bonds owned by old owned atoms

  int *maxhalf;                 // bond index for maxstrain bond of old atoms

  int *eligible;                // 0/1 flag for bias on one of old atom's bonds

  double *maxhalfstrain;        // strain value for maxstrain bond of old atoms

  struct OneBond {             // single IJ bond, atom I is owner

    double r0;                 // original relaxed bond length

    double boostcoeff;         // boost coefficient

    tagint jtag;               // global index of J atom in bond IJ

    int j;                     // local index of J atom in bond IJ

  };

  int *old2now;                 // o2n[i] = current local index of old atom I

                                // may be -1 if ghost atom has drifted

  tagint *tagold;               // IDs of atoms when bonds were formed

                                // 0 if a ghost atom is not in Dcut neigh list

  double **xold;                // coords of atoms when bonds were formed

  // vectors used to find maxstrain bonds within a local domain

  struct OneBond **bonds;      // 2d array of bonds for owned atoms

  int *numbond;                // number of bonds for each owned atom

  int maxatom;                 // size of these vectors, nlocal + nghost

  double *maxstrain;           // max-strain of any bond atom I is part of

                               //   for owned and ghost atoms

  double *maxstrain_region;    // max-strain of any neighbor atom J of atom I

  double *maxstrain_domain;    // max-strain of any neighbor atom J of atom I

                               //   for owned and ghost atoms

  int *maxstrain_bondindex;    // index of max-strain bond of each atom I

                               //   just for owned atoms

  tagint *biasflag;            // atoms in biased bonds marked with bond partner

                               //   for owned and ghost atoms

  // list of boosted bonds that this proc will bias

  // data struct used to persist biascoeffs when bond list is re-created

  int maxboost;                // allocated size of boost list

  int nboost;                  // # of boosted bonds I own

  int *boost;                  // index of atom I in each boosted bond

  struct OneCoeff {

    double biascoeff;

    tagint jtag;

  };

  struct OneCoeff **clist;     // list of bond coeffs for each atom's bonds

  int *numcoeff;               // # of coeffs per atom

  int maxcoeff;                // allocate size of clist

  int maxcoeffperatom;         // allocated # of columns in clist

  // list of biased bonds this proc owns

  int maxbias;                 // allocated size of bias list

  int nbias;                   // # of biased bonds I own

  int *bias;                   // index of biased bonds in my bond list

  // extra timers

  //double timefirst,timesecond,timethird,timefourth;

  //double timefifth,timesixth,timeseventh,timetotal;

  // histogramming of bond boost cooeficients

  // private methods

  int histo_every,histo_count,histo_print,histo_steps;

  double histo_delta,invhisto_delta,histo_lo;

  bigint *histo,*allhisto;

  void grow_bond();

  void grow_coeff();

};

}