From eea30c5b766f15f2c011ba0736c6bf4b4c46dfc1 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 26 Feb 2019 16:17:31 -0700 Subject: [PATCH] half-bond lists and no bond migration for local hyper --- src/REPLICA/fix_hyper_global.cpp | 124 ++- src/REPLICA/fix_hyper_global.h | 4 +- src/REPLICA/fix_hyper_local.cpp | 1299 +++++++++++++++--------------- src/REPLICA/fix_hyper_local.h | 142 ++-- 4 files changed, 814 insertions(+), 755 deletions(-) diff --git a/src/REPLICA/fix_hyper_global.cpp b/src/REPLICA/fix_hyper_global.cpp index d235b06dc0..e43f1431a9 100644 --- a/src/REPLICA/fix_hyper_global.cpp +++ b/src/REPLICA/fix_hyper_global.cpp @@ -32,12 +32,12 @@ 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() /* ---------------------------------------------------------------------- */ @@ -102,7 +102,6 @@ int FixHyperGlobal::setmask() { int mask = 0; mask |= PRE_NEIGHBOR; - mask |= PRE_FORCE; mask |= PRE_REVERSE; mask |= THERMO_ENERGY; return mask; @@ -124,6 +123,10 @@ void FixHyperGlobal::init() 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 @@ -167,15 +170,56 @@ void FixHyperGlobal::setup_pre_reverse(int eflag, int vflag) 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++) { @@ -196,6 +240,8 @@ void FixHyperGlobal::pre_neighbor() } if (flag) error->one(FLERR,"Fix hyper/global bond atom not found"); + + */ } /* ---------------------------------------------------------------------- */ @@ -204,15 +250,16 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) { 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; @@ -225,8 +272,8 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) 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; @@ -238,7 +285,7 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) // 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; @@ -255,25 +302,34 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) 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; @@ -281,13 +337,14 @@ void FixHyperGlobal::pre_reverse(int /* eflag */, int /* vflag */) 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]; @@ -364,13 +421,16 @@ void FixHyperGlobal::build_bond_list(int natom) 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]; @@ -386,14 +446,11 @@ void FixHyperGlobal::build_bond_list(int natom) 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"); } /* ---------------------------------------------------------------------- */ @@ -419,7 +476,7 @@ double FixHyperGlobal::compute_vector(int i) // 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 @@ -438,8 +495,9 @@ double FixHyperGlobal::compute_vector(int i) 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; } diff --git a/src/REPLICA/fix_hyper_global.h b/src/REPLICA/fix_hyper_global.h index d8b1cef425..42dd64e145 100644 --- a/src/REPLICA/fix_hyper_global.h +++ b/src/REPLICA/fix_hyper_global.h @@ -76,10 +76,12 @@ class FixHyperGlobal : public FixHyper { // 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 diff --git a/src/REPLICA/fix_hyper_local.cpp b/src/REPLICA/fix_hyper_local.cpp index 60198f4ddf..99dd1945ad 100644 --- a/src/REPLICA/fix_hyper_local.cpp +++ b/src/REPLICA/fix_hyper_local.cpp @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- + /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -33,33 +33,30 @@ using namespace LAMMPS_NS; using namespace FixConst; -#define DELTABOOST 16 -#define BOOSTINIT 1.0 +#define DELTABOND 16384 +#define DELTABIAS 16 +#define COEFFINIT 1.0 #define COEFFMAX 1.2 +#define MAXBONDPERATOM 30 #define BIG 1.0e20 -enum{STRAIN,STRAINREGION,BIASFLAG}; +enum{STRAIN,STRAINDOMAIN,BIASFLAG,BIASCOEFF}; enum{IGNORE,WARN,ERROR}; /* ---------------------------------------------------------------------- */ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : FixHyper(lmp, narg, arg), old2now(NULL), xold(NULL), tagold(NULL), - bonds(NULL), numbond(NULL), maxstrain(NULL), maxstrain_region(NULL), - maxstrain_bondindex(NULL), biasflag(NULL), boost(NULL), - histo(NULL), allhisto(NULL) + blist(NULL), maxstrain(NULL), maxstrain_domain(NULL), + biasflag(NULL), bias(NULL) { + // NOTE: need to add vecs/arrays to constructor list + // error checks - // solution for tagint != int is to worry about storing - // local index vs global ID in same variable - // maybe need to declare them all tagint, not int if (atom->map_style == 0) error->all(FLERR,"Fix hyper/local command requires atom map"); - if (sizeof(tagint) != sizeof(int)) - error->all(FLERR,"Fix hyper/local requires tagint = int"); - // parse args if (narg < 10) error->all(FLERR,"Illegal fix hyper/local command"); @@ -79,10 +76,10 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : tequil = force->numeric(FLERR,arg[6]); dcut = force->numeric(FLERR,arg[7]); alpha_user = force->numeric(FLERR,arg[8]); - boosttarget = force->numeric(FLERR,arg[9]); + boost_target = force->numeric(FLERR,arg[9]); if (cutbond < 0.0 || qfactor < 0.0 || vmax < 0.0 || - tequil <= 0.0 || dcut <= 0.0 || alpha_user <= 0.0 || boosttarget < 1.0) + tequil <= 0.0 || dcut <= 0.0 || alpha_user <= 0.0 || boost_target < 1.0) error->all(FLERR,"Illegal fix hyper/local command"); invqfactorsq = 1.0 / (qfactor*qfactor); @@ -92,36 +89,11 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : // optional args - histoflag = 0; - lostbond = IGNORE; checkbias = 0; - checkcoeff = 0; int iarg = 10; while (iarg < narg) { - /* NOTE: do not enable this yet, need to think about it differently - if (strcmp(arg[iarg],"histo") == 0) { - if (iarg+5 > narg) error->all(FLERR,"Illegal fix hyper/local command"); - histoflag = 1; - histo_every = force->inumeric(FLERR,arg[iarg+1]); - histo_count = force->inumeric(FLERR,arg[iarg+2]); - histo_delta = force->numeric(FLERR,arg[iarg+3]); - histo_print = force->inumeric(FLERR,arg[iarg+4]); - if (histo_every <= 0 || histo_count % 2 || - histo_delta <= 0.0 || histo_print <= 0) - error->all(FLERR,"Illegal fix hyper/local command"); - iarg += 5; - */ - - if (strcmp(arg[iarg],"lostbond") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix hyper/local command"); - if (strcmp(arg[iarg+1],"error") == 0) lostbond = ERROR; - else if (strcmp(arg[iarg+1],"warn") == 0) lostbond = WARN; - else if (strcmp(arg[iarg+1],"ignore") == 0) lostbond = IGNORE; - else error->all(FLERR,"Illegal fix hyper/local command"); - iarg += 2; - - } else if (strcmp(arg[iarg],"check/bias") == 0) { + if (strcmp(arg[iarg],"check/bias") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal fix hyper/local command"); checkbias = 1; checkbias_every = force->inumeric(FLERR,arg[iarg+1]); @@ -131,56 +103,51 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : else error->all(FLERR,"Illegal fix hyper/local command"); iarg += 3; - } else if (strcmp(arg[iarg],"check/coeff") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix hyper/local command"); - checkcoeff = 1; - checkcoeff_every = force->inumeric(FLERR,arg[iarg+1]); - if (strcmp(arg[iarg+2],"error") == 0) checkcoeff_flag = ERROR; - else if (strcmp(arg[iarg+2],"warn") == 0) checkcoeff_flag = WARN; - else if (strcmp(arg[iarg+2],"ignore") == 0) checkcoeff_flag = IGNORE; - else error->all(FLERR,"Illegal fix hyper/local command"); - iarg += 3; - } else error->all(FLERR,"Illegal fix hyper/local command"); } // per-atom data structs maxbond = 0; - bonds = NULL; - numbond = NULL; + blist = NULL; + + maxatom = 0; maxstrain = NULL; - maxstrain_region = NULL; - maxstrain_bondindex = NULL; + maxstrain_domain = NULL; biasflag = NULL; - maxold = old_nall = 0; - old2now = NULL; - xold = NULL; // NOTE: don't really need except to monitor drift - tagold = NULL; + maxlocal = nlocal_old = 0; + numbond = NULL; + maxhalf = NULL; + eligible = NULL; + maxhalfstrain = NULL; - nboost = maxboost = 0; - boost = NULL; + maxall = nall_old = 0; + xold = NULL; + tagold = NULL; + old2now = NULL; + + nbias = maxbias = 0; + bias = NULL; + + maxcoeff = 0; + maxcoeffperatom = 0; + numcoeff = NULL; + clist = NULL; // maxbondperatom = max # of bonds any atom is part of // will be reset in bond_build() - // set comm size needed by this fix + // set comm sizes needed by this fix + // NOTE: remove MBPA when minimize reverse Cij comm maxbondperatom = 1; comm_forward = 1; - comm_reverse = 1; - - // perform initial allocation of atom-based arrays - // register with Atom class - - grow_arrays(atom->nmax); - atom->add_callback(0); + comm_reverse = MAXBONDPERATOM; me = comm->me; firstflag = 1; - allbonds = 0; - allboost = 0.0; + allbias = 0.0; starttime = update->ntimestep; nostrainyet = 1; @@ -188,36 +155,31 @@ FixHyperLocal::FixHyperLocal(LAMMPS *lmp, int narg, char **arg) : nevent = 0; nevent_atom = 0; mybias = 0.0; - - histo = allhisto = NULL; - if (histoflag) { - invhisto_delta = 1.0 / histo_delta; - histo_lo = 1.0 - (histo_count/2 * histo_delta); - histo = new bigint[histo_count+2]; - allhisto = new bigint[histo_count+2]; - } } /* ---------------------------------------------------------------------- */ FixHyperLocal::~FixHyperLocal() { - memory->destroy(bonds); - memory->destroy(numbond); - - atom->delete_callback(id,0); + memory->destroy(blist); memory->destroy(maxstrain); - memory->destroy(maxstrain_region); - memory->destroy(maxstrain_bondindex); + memory->destroy(maxstrain_domain); memory->destroy(biasflag); - memory->destroy(old2now); + memory->destroy(numbond); + memory->destroy(maxhalf); + memory->destroy(eligible); + memory->destroy(maxhalfstrain); + memory->destroy(xold); memory->destroy(tagold); - memory->destroy(boost); - delete [] histo; - delete [] allhisto; + memory->destroy(old2now); + + memory->destroy(bias); + + memory->destroy(numcoeff); + memory->destroy(clist); } /* ---------------------------------------------------------------------- */ @@ -237,24 +199,19 @@ int FixHyperLocal::setmask() void FixHyperLocal::init_hyper() { ghost_toofar = 0; - lostbond_partner = 0; - lostbond_coeff = 0.0; checkbias_count = 0; - checkcoeff_count = 0; maxdriftsq = 0.0; maxbondlen = 0.0; - maxboostcoeff = 0.0; - minboostcoeff = BIG; - sumboostcoeff = 0.0; - nboost_running = 0; + maxbiascoeff = 0.0; + minbiascoeff = BIG; + sumbiascoeff = 0.0; + nbias_running = 0; nobias_running = 0; rmaxever = 0.0; rmaxeverbig = 0.0; nbondbuild = 0; time_bondbuild = 0.0; - - if (histoflag) histo_steps = 0; } /* ---------------------------------------------------------------------- */ @@ -262,7 +219,7 @@ void FixHyperLocal::init_hyper() void FixHyperLocal::init() { // for newton off, bond force bias will not be applied correctly - // bonds that straddle 2 procs + // for bonds that straddle 2 procs // warn if molecular system, since near-neighbors may not appear in neigh list // user should not be including bonded atoms as hyper "bonds" @@ -271,7 +228,7 @@ void FixHyperLocal::init() if (atom->molecular && me == 0) error->warning(FLERR,"Hyper local for molecular systems " - "requires care in defining hyperdynamics bonds"); + "requires care in defining hyperdynamic bonds"); // cutghost = communication cutoff as calculated by Neighbor and Comm // error if cutghost is smaller than Dcut @@ -285,30 +242,41 @@ void FixHyperLocal::init() cutghost = comm->cutghostuser; if (cutghost < dcut) - error->all(FLERR,"Fix hyper/local bond cutoff exceeds ghost atom range - " + error->all(FLERR,"Fix hyper/local domain cutoff exceeds ghost atom range - " "use comm_modify cutoff command"); if (cutghost < dcut+cutbond/2.0 && me == 0) error->warning(FLERR,"Fix hyper/local ghost atom range " "may not allow for atom drift between events"); } - alpha = update->dt / alpha_user; - // need an occasional full neighbor list with cutoff = Dcut + // need occasional full neighbor list with cutoff = Dcut + // used for finding maxstrain of neighbor bonds out to Dcut // do not need to include neigh skin in cutoff, // b/c this list will be built every time bond_build() is called // NOTE: what if pair style list cutoff > Dcut // or what if neigh skin is huge? - 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; - neighbor->requests[irequest]->cut = 1; - neighbor->requests[irequest]->cutoff = dcut; - neighbor->requests[irequest]->occasional = 1; + int irequest_full = neighbor->request(this,instance_me); + neighbor->requests[irequest_full]->pair = 0; + neighbor->requests[irequest_full]->fix = 1; + neighbor->requests[irequest_full]->half = 0; + neighbor->requests[irequest_full]->full = 1; + neighbor->requests[irequest_full]->cut = 1; + neighbor->requests[irequest_full]->cutoff = dcut; + neighbor->requests[irequest_full]->occasional = 1; + + // also need occasional half neighbor list derived from pair style + // used for building local bond list + // no specified cutoff, should be longer than cutbond + // this list will also be built (or derived/copied) + // every time bond_build() is called + + int irequest_half = neighbor->request(this,instance_me); + neighbor->requests[irequest_half]->pair = 0; + neighbor->requests[irequest_half]->fix = 1; + neighbor->requests[irequest_half]->occasional = 1; // extra timing output @@ -318,16 +286,18 @@ void FixHyperLocal::init() /* ---------------------------------------------------------------------- */ -void FixHyperLocal::init_list(int /* id */, NeighList *ptr) +void FixHyperLocal::init_list(int id, NeighList *ptr) { - list = ptr; + if (id == 1) listfull = ptr; + else if (id == 2) listhalf = ptr; } /* ---------------------------------------------------------------------- */ void FixHyperLocal::setup_pre_neighbor() { - // called for dynamics and minimization NOTE: check if min is needed? + // called for dynamics and minimization + // NOTE: check if needed for min, I think so b/c of Cij persist? pre_neighbor(); } @@ -337,8 +307,8 @@ void FixHyperLocal::setup_pre_neighbor() void FixHyperLocal::setup_pre_reverse(int eflag, int vflag) { // only called for dynamics, not minimization - // setupflag prevents boostostat update of boost coeffs in setup - // also prevents increments of nboost_running, nbias_running, sumboostcoeff + // setupflag prevents boostostat update of bias coeffs in setup + // also prevents increments of nbias_running, nobias_running, sumbiascoeff setupflag = 1; pre_reverse(eflag,vflag); @@ -349,67 +319,67 @@ void FixHyperLocal::setup_pre_reverse(int eflag, int vflag) void FixHyperLocal::pre_neighbor() { - int i,m,n,ilocal,jlocal; + int i,m,iold,jold,ilocal,jlocal; + double distsq; - // convert global ID bond partners back to local indices - // need to use closest_image() so can calculate bond lengths - // error flag should not happen for a well-behaved system - // b/c are only looking up bond partners inside or near my sub-domain + // 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; - int nlocal = atom->nlocal; - int missing = 0; - double missing_coeff = 0.0; + for (m = 0; m < nblocal; m++) { + iold = blist[m].iold; + jold = blist[m].jold; + ilocal = old2now[iold]; + jlocal = old2now[jold]; - for (i = 0; i < nlocal; i++) { - n = numbond[i]; - for (m = 0; m < n; m++) { - jlocal = atom->map(bonds[i][m].jtag); - if (jlocal >= 0) bonds[i][m].j = domain->closest_image(i,jlocal); - else { - bonds[i][m].j = -1; - missing++; - missing_coeff += bonds[i][m].boostcoeff; - if (lostbond != IGNORE) { - char str[128]; - sprintf(str,"Fix hyper/local bond info missing for bond " - TAGINT_FORMAT "," TAGINT_FORMAT - " with coeff %g at step " BIGINT_FORMAT, - atom->tag[i],bonds[i][m].jtag,bonds[i][m].boostcoeff, - update->ntimestep); - if (lostbond == ERROR) error->one(FLERR,str); - if (lostbond == WARN) error->warning(FLERR,str); - } - } + if (ilocal < 0) { + ilocal = atom->map(tagold[iold]); + ilocal = domain->closest_image(xold[iold],ilocal); + if (ilocal < 0) + error->one(FLERR,"Fix hyper/local 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/local 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; } - lostbond_partner += missing; - lostbond_coeff += missing_coeff; - - // set old2now to point to current local atom indices + // set remaining old2now values to point to current local atom indices + // if old2now >= 0, already set by bond loop above // only necessary for tagold entries > 0 // because if tagold = 0, atom is not active in Dcut neighbor list // must be done after atoms migrate and ghost atoms setup via comm->borders() - // does not matter if there are multiple ghost copies of a global ID - // since only need the ghost atom strain, not its coordinates - // NOTE: maybe need not use closest image, b/c old2now only used to access - // maxstrain values, which will be same for any copy of ghost atom ?? - // also need it to compute maxdriftsq correctly when a proc spans a PBC + // does not matter which atom (owned or ghost) that atom->map() finds + // b/c old2now is only used to access maxstrain() or biasflag() + // which will be identical for every copy of the same atom ID - double distsq; - - for (i = 0; i < old_nall; i++) { + for (i = 0; i < nall_old; i++) { + if (old2now[i] >= 0) continue; if (tagold[i] == 0) continue; ilocal = atom->map(tagold[i]); - ilocal = domain->closest_image(xold[i],ilocal); old2now[i] = ilocal; - - if (ilocal >= 0) { - distsq = MathExtra::distsq3(x[ilocal],xold[i]); - maxdriftsq = MAX(distsq,maxdriftsq); - } else ghost_toofar++; + if (ilocal < 0) ghost_toofar++; } } @@ -417,93 +387,116 @@ void FixHyperLocal::pre_neighbor() void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) { - int i,j,m,ii,jj,inum,jnum,iold,jold,nbond,bondindex; - tagint itag,jtag; + int i,j,m,ii,jj,inum,jnum,iold,jold,ibond,nbond,ijhalf,ncount; double xtmp,ytmp,ztmp,delx,dely,delz; - double r,r0,estrain,emax,vbias,fbias,fbiasr,boostcoeff; + double r,r0,estrain,emax,ebias,vbias,fbias,fbiasr,biascoeff; + double halfstrain,selfstrain; int *ilist,*jlist,*numneigh,**firstneigh; //double time1,time2,time3,time4,time5,time6,time7,time8; //time1 = MPI_Wtime(); - // compute current maxstrain and maxstrain_bond for each owned atom - // use per-atom full bond list - // this is double-calculating for IJ and JI bonds - // could compute once, but would have to find/store index of JI bond - // order two I,J atoms consistently for IJ and JI calcs - // to insure no round-off issue when comparing maxstrain values of I,J + // reallocate local vectors if necessary + + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + if (maxatom < nall) { + maxatom = atom->nmax; + memory->grow(maxstrain,maxatom,"hyper/local:maxstrain"); + memory->grow(maxstrain_domain,maxatom,"hyper/local:maxstrain_domain"); + if (checkbias) memory->grow(biasflag,maxatom,"hyper/local:biasflag"); + } + + // each old atom I's owned bond with max strain is eligible for biasing + + for (iold = 0; iold < nlocal_old; iold++) eligible[iold] = 1; + + // ------------------------------------------------------------- + // stage 1: + // maxstrain[i] = max abs value of strain of any bond atom I is part of + // reverse/forward comm so know it for all current owned and ghost atoms + // ------------------------------------------------------------- + + // compute estrain = current abs value strain of each owned bond + // blist = bondlist from last event + // mark atom I ineligible if it has no bonds + // also store: + // maxhalf = which owned bond is maxstrain for each old atom I + // maxhalfstrain = strain of that bond for each old atom I + + for (i = 0; i < nall; i++) maxstrain[i] = 0.0; double **x = atom->x; - tagint *tag = atom->tag; - int nlocal = atom->nlocal; - int mybonds = 0; - nostrainyet = 0; - - for (i = 0; i < nlocal; i++) { - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itag = tag[i]; - emax = 0.0; - bondindex = -1; - jtag = 0; - - nbond = numbond[i]; - mybonds += nbond; - for (m = 0; m < nbond; m++) { - j = bonds[i][m].j; - if (j < 0) continue; - jtag = bonds[i][m].jtag; - r0 = bonds[i][m].r0; - if (itag < jtag) { - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - } else { - delx = x[j][0] - xtmp; - dely = x[j][1] - ytmp; - delz = x[j][2] - ztmp;; - } + m = 0; + for (iold = 0; iold < nlocal_old; iold++) { + nbond = numbond[iold]; + if (!nbond) { + eligible[iold] = 0; + continue; + } + halfstrain = 0.0; + for (ibond = 0; ibond < nbond; ibond++) { + i = blist[m].i; + j = blist[m].j; + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; r = sqrt(delx*delx + dely*dely + delz*delz); maxbondlen = MAX(r,maxbondlen); + r0 = blist[m].r0; estrain = fabs(r-r0) / r0; - if (estrain > emax) { - emax = estrain; - bondindex = m; + maxstrain[i] = MAX(maxstrain[i],estrain); + maxstrain[j] = MAX(maxstrain[j],estrain); + if (estrain > halfstrain) { + halfstrain = estrain; + ijhalf = m; } + m++; } - maxstrain[i] = emax; - maxstrain_bondindex[i] = bondindex; + maxhalf[iold] = ijhalf; + maxhalfstrain[iold] = halfstrain; } //time2 = MPI_Wtime(); - // forward comm to acquire maxstrain of all ghost atoms + // reverse comm acquires maxstrain of all current owned atoms + // needed b/c only saw half the bonds of each atom + // also needed b/c bond list may refer to old owned atoms that are now ghost + // forward comm acquires maxstrain of all current ghost atoms commflag = STRAIN; + comm->reverse_comm_fix(this); comm->forward_comm_fix(this); //time3 = MPI_Wtime(); - // use original Dcut neighbor list to check maxstrain of all neighbor atoms - // set maxstrain_region of I atoms = maxstrain of I and all J neighs - // neighbor list has old indices for IJ b/c reneighboring may have occurred + // ------------------------------------------------------------- + // stage 2: + // maxstrain_domain[i] = maxstrain of atom I and all its J neighs out to Dcut + // reverse/forward comm so know it for all current owned and ghost atoms + // ------------------------------------------------------------- + + // use full Dcut neighbor list to check maxstrain of all neighbor atoms + // NOTE: is II loop the same as iold over nlocal_old ?? + // neighlist is from last event + // has old indices for I,J (reneighboring may have occurred) // use old2now[] to convert to current indices - // if neighbor is not currently known (too far away), - // then assume it was part of an event and its strain = qfactor - // this double loop sets maxstrain_region of mostly owned atoms - // but possibly some ghost atoms as well + // if J is unknown (drifted ghost), + // assume it was part of an event and its strain = qfactor + // mark atom I ineligible for biasing if: + // its maxhalfstrain < maxstrain (J atom owns the IJ bond) + // its maxstrain < maxstrain_domain + // ncount > 1 (break tie by making all atoms with tie value ineligible) + // if ncount > 1, also flip sign of maxstrain_domain for atom I - int nall = nlocal + atom->nghost; - for (i = 0; i < nall; i++) maxstrain_region[i] = 0.0; + for (i = 0; i < nall; i++) maxstrain_domain[i] = 0.0; - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // find largest distance from subbox that a ghost atom is with strain < qfactor + inum = listfull->inum; + ilist = listfull->ilist; + numneigh = listfull->numneigh; + firstneigh = listfull->firstneigh; double rmax = rmaxever; double rmaxbig = rmaxeverbig; @@ -515,23 +508,36 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) jlist = firstneigh[iold]; jnum = numneigh[iold]; - // I and J may be ghost atoms + // I or J may be ghost atoms // will always know I b/c atoms do not drift that far // but may no longer know J if hops outside cutghost // in that case, assume it performed an event, its strain = qfactor + // this assumes cutghost is sufficiently longer than Dcut i = old2now[iold]; - emax = maxstrain[i]; - + emax = selfstrain = maxstrain[i]; + ncount = 0; + for (jj = 0; jj < jnum; jj++) { jold = jlist[jj]; j = old2now[jold]; - if (j >= 0) emax = MAX(emax,maxstrain[j]); - else { + + // special case for missing (drifted) J atom + + if (j < 0) { emax = MAX(emax,qfactor); + if (selfstrain == qfactor) ncount++; continue; } + emax = MAX(emax,maxstrain[j]); + if (selfstrain == maxstrain[j]) ncount++; + + // diagnostic + // tally largest distance from subbox that a ghost atom is (rmaxbig) + // and the largest distance if strain < qfactor (rmax) + // NOTE: could this be removed from loop ?? + if (j >= nlocal) { if (x[j][0] < sublo[0]) rmaxbig = MAX(rmaxbig,sublo[0]-x[j][0]); if (x[j][1] < sublo[1]) rmaxbig = MAX(rmaxbig,sublo[1]-x[j][1]); @@ -550,9 +556,17 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) } } - maxstrain_region[i] = emax; + if (maxhalfstrain[iold] < selfstrain) eligible[iold] = 0; + if (selfstrain < emax) eligible[iold] = 0; + else if (ncount > 1) { + eligible[iold] = 0; + emax = -emax; + } + maxstrain_domain[i] = emax; } + // diagnostic // NOTE: optional, should skip + double rmax2[2],rmax2all[2]; rmax2[0] = rmax; rmax2[1] = rmaxbig; @@ -560,80 +574,76 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) rmaxever = rmax2all[0]; rmaxeverbig = rmax2all[1]; - MPI_Allreduce(&mybonds,&allbonds,1,MPI_INT,MPI_SUM,world); - //time4 = MPI_Wtime(); - // reverse comm to acquire maxstrain_region from ghost atoms - // needed b/c neighbor list referred to old owned atoms, - // so above loop may set maxstrain_region of ghost atoms - // forward comm to acquire maxstrain_region of all ghost atoms + // reverse comm to acquire maxstrain_domain from ghost atoms + // needed b/c neigh list may refer to old owned atoms that are now ghost + // forward comm acquires maxstrain_domain of all current ghost atoms - commflag = STRAINREGION; + commflag = STRAINDOMAIN; comm->reverse_comm_fix(this); comm->forward_comm_fix(this); //time5 = MPI_Wtime(); - // identify biased bonds and add to boost list - // for each max-strain bond IJ of atom I: - // bias this bond only if all these conditions hold: - // itag < jtag, so bond is only biased once - // maxstrain[i] = maxstrain_region[i] - // maxstrain[j] = maxstrain_region[j] - // NOTE: also need to check that maxstrain[i] = maxstrain[j] ?? - // don't think so, b/c maxstrain_region[i] includes maxstrain[i] + // ------------------------------------------------------------- + // stage 3: + // create bias = list of Nbias biased bonds this proc owns + // ------------------------------------------------------------- - nboost = 0; + // identify biased bonds and add to bias list + // bias the I,J maxhalf bond of atom I only if all these conditions hold: + // maxstrain[i] = maxstrain_domain[i] (checked in stage 2) + // maxstrain[j] = maxstrain_domain[j] (checked here) + // I is not part of an I,J bond with > strain owned by some J (checked in 2) + // no ties with other maxstrain bonds in atom I's domain (chedcked in 2) - for (i = 0; i < nlocal; i++) { - if (numbond[i] == 0) continue; - itag = tag[i]; - j = bonds[i][maxstrain_bondindex[i]].j; - jtag = tag[j]; - if (itag > jtag) continue; - - if (maxstrain[i] != maxstrain_region[i]) continue; - if (maxstrain[j] != maxstrain_region[j]) continue; - - if (nboost == maxboost) { - maxboost += DELTABOOST; - memory->grow(boost,maxboost,"hyper/local:boost"); + nbias = 0; + for (iold = 0; iold < nlocal_old; iold++) { + if (eligible[iold] == 0) continue; + j = blist[maxhalf[iold]].j; + if (maxstrain[j] != maxstrain_domain[j]) continue; + if (nbias == maxbias) { + maxbias += DELTABIAS; + memory->grow(bias,maxbias,"hyper/local:bias"); } - boost[nboost++] = i; + bias[nbias++] = ibond; } //time6 = MPI_Wtime(); - // apply boost force to bonds with locally max strain + // ------------------------------------------------------------- + // stage 4: + // apply bias force to bonds with locally max strain + // ------------------------------------------------------------- double **f = atom->f; int nobias = 0; mybias = 0.0; - for (int iboost = 0; iboost < nboost; iboost++) { - i = boost[iboost]; - emax = maxstrain[i]; - if (emax >= qfactor) { + for (int ibias = 0; ibias < nbias; ibias++) { + m = bias[ibias]; + i = blist[m].i; + j = blist[m].j; + + if (maxstrain[i] >= qfactor) { nobias++; continue; } - m = maxstrain_bondindex[i]; - j = bonds[i][m].j; - r0 = bonds[i][m].r0; - boostcoeff = bonds[i][m].boostcoeff; - - vbias = boostcoeff * vmax * (1.0 - emax*emax*invqfactorsq); - fbias = boostcoeff * 2.0 * vmax * emax / (qfactor*qfactor * r0); - delx = x[i][0] - x[j][0]; dely = x[i][1] - x[j][1]; delz = x[i][2] - x[j][2]; r = sqrt(delx*delx + dely*dely + delz*delz); - fbiasr = fbias / r; + r0 = blist[m].r0; + //ebias = (r-r0) / r0; + ebias = fabs(r-r0) / r0; + biascoeff = blist[m].biascoeff; + vbias = biascoeff * vmax * (1.0 - ebias*ebias*invqfactorsq); + fbias = biascoeff * 2.0 * vmax * ebias * invqfactorsq; + fbiasr = fbias / r0 / r; f[i][0] += delx*fbiasr; f[i][1] += dely*fbiasr; f[i][2] += delz*fbiasr; @@ -647,111 +657,72 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) //time7 = MPI_Wtime(); - // no boostostat update of boost coeffs when pre_reverse called from setup() - // nboost_running, nobias_running, sumboostcoeff only incremented on run steps + // ------------------------------------------------------------- + // apply boostostat to bias coeffs of all bonds I own + // ------------------------------------------------------------- + + // no boostostat update when pre_reverse called from setup() + // nbias_running, nobias_running, sumbiascoeff only incremented on run steps // NOTE: maybe should also not bias any bonds on firststep of this fix if (setupflag) return; - nboost_running += nboost; + nbias_running += nbias; nobias_running += nobias; - // apply boostostat to boost coefficients of all bonds of all owned atoms - // use per-atom full bond list - // this is double-calculating for IJ and JI bonds - // should be identical for both, b/c emax is the same - // could compute once, but would have to find/store index of JI bond - // delta in boost coeff is function of maxboost_region vs target boost - // maxboost_region is function of two maxstrain_regions for I,J - // NOTE: if J is lost to I but not vice versa, then biascoeff IJ != JI + // loop over bonds I own to adjust bias coeff + // delta in boost coeff is function of maxboost_domain vs target boost + // maxboost_domain is function of two maxstrain_domains for I,J - double myboost = 0.0; - double emaxi,emaxj,maxboost_region; + double mybias = 0.0; + double emaxi,emaxj,maxboost_domain; - for (i = 0; i < nlocal; i++) { - emaxi = maxstrain_region[i]; - nbond = numbond[i]; - for (m = 0; m < nbond; m++) { - j = bonds[i][m].j; - if (j < 0) continue; - emaxj = maxstrain_region[j]; - emax = MAX(emaxi,emaxj); - if (emax < qfactor) vbias = vmax * (1.0 - emax*emax*invqfactorsq); - else vbias = 0.0; - boostcoeff = bonds[i][m].boostcoeff; - maxboost_region = exp(beta * boostcoeff*vbias); - boostcoeff -= alpha * (maxboost_region-boosttarget) / boosttarget; - // COMMENT OUT for now - need better way to bound boostcoeff - //boostcoeff = MIN(boostcoeff,COEFFMAX); - myboost += boostcoeff; - maxboostcoeff = MAX(maxboostcoeff,boostcoeff); - minboostcoeff = MIN(minboostcoeff,boostcoeff); - bonds[i][m].boostcoeff = boostcoeff; - } + for (m = 0; m < nblocal; m++) { + i = blist[m].i; + j = blist[m].j; + emaxi = fabs(maxstrain_domain[i]); + emaxj = fabs(maxstrain_domain[j]); + emax = MAX(emaxi,emaxj); + if (emax < qfactor) vbias = vmax * (1.0 - emax*emax*invqfactorsq); + else vbias = 0.0; + + biascoeff = blist[m].biascoeff; + maxboost_domain = exp(beta * biascoeff*vbias); + biascoeff -= alpha * (maxboost_domain-boost_target) / boost_target; + blist[m].biascoeff = biascoeff; + + // stats + + mybias += biascoeff; + maxbiascoeff = MAX(maxbiascoeff,biascoeff); + minbiascoeff = MIN(minbiascoeff,biascoeff); } // running stats - MPI_Allreduce(&myboost,&allboost,1,MPI_DOUBLE,MPI_SUM,world); - if (allbonds) sumboostcoeff += allboost/allbonds; + MPI_Allreduce(&mybias,&allbias,1,MPI_DOUBLE,MPI_SUM,world); + if (allbonds) sumbiascoeff += allbias/allbonds; - // histogram the bond coeffs and output as requested - // do not double count each bond + // ------------------------------------------------------------- + // extra diagnostics if requested + // ------------------------------------------------------------- - if (histoflag && update->ntimestep % histo_every == 0) { - if (histo_steps == 0) - for (i = 0; i < histo_count+2; i++) histo[i] = 0; - histo_steps++; - - int ihisto; - for (i = 0; i < nlocal; i++) { - nbond = numbond[i]; - for (m = 0; m < nbond; m++) { - if (tag[i] > bonds[i][m].jtag) continue; - boostcoeff = bonds[i][m].boostcoeff; - if (boostcoeff < histo_lo) ihisto = -1; - else ihisto = static_cast ((boostcoeff-histo_lo) * invhisto_delta); - if (ihisto >= histo_count) ihisto = histo_count; - histo[ihisto+1]++; - } - } - - if (update->ntimestep % histo_print == 0) { - MPI_Allreduce(histo,allhisto,histo_count+2,MPI_LMP_BIGINT,MPI_SUM,world); - - bigint total = 0; - for (i = 0; i < histo_count+2; i++) total += allhisto[i]; - - if (me == 0) { - if (screen) { - fprintf(screen,"Histogram of bias coeffs:\n"); - for (i = 0; i < histo_count+2; i++) - fprintf(screen," %g",1.0*allhisto[i]/total); - fprintf(screen,"\n"); - } - if (logfile) { - fprintf(logfile,"Histogram of bias coeffs:\n"); - for (i = 0; i < histo_count+2; i++) - fprintf(logfile," %g",1.0*allhisto[i]/total); - fprintf(logfile,"\n"); - } - } - } - } - - // check for any biased bonds that are too close to each other + // if requsted, check for any biased bonds that are too close to each other // keep a running count for output + // requires 2 additional local comm operations if (checkbias && update->ntimestep % checkbias_every == 0) { // mark each atom in a biased bond with ID of partner - // nboost loop will mark some ghost atoms + // nbias loop will mark some ghost atoms for (i = 0; i < nall; i++) biasflag[i] = 0; - for (int iboost = 0; iboost < nboost; iboost++) { - i = boost[iboost]; - m = maxstrain_bondindex[i]; - j = bonds[i][m].j; + tagint *tag = atom->tag; + + for (int ibias = 0; ibias < nbias; ibias++) { + m = bias[ibias]; + i = blist[m].i; + j = blist[m].j; biasflag[i] = tag[j]; biasflag[j] = tag[i]; } @@ -766,7 +737,7 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) // I and J may be ghost atoms // only continue if I is a biased atom - // if J is unknonw (drifted ghost) just ignore + // if J is unknown (drifted ghost) just ignore // if J is biased and is not bonded to I, then flag as too close for (ii = 0; ii < inum; ii++) { @@ -785,32 +756,6 @@ void FixHyperLocal::pre_reverse(int /* eflag */, int /* vflag */) } } } - - // check for any bond bias coeffcients that do not match - // cannot check unless both atoms IJ are owned by this proc - // keep a running count for output - - if (checkcoeff && update->ntimestep % checkcoeff_every == 0) { - int jb,jbonds; - - for (i = 0; i < nlocal; i++) { - nbond = numbond[i]; - for (m = 0; m < nbond; m++) { - if (tag[i] > bonds[i][m].jtag) continue; - j = bonds[i][m].j; - if (j < 0) continue; - if (j >= nlocal) continue; - itag = tag[i]; - jbonds = numbond[j]; - for (jb = 0; jb < jbonds; jb++) - if (bonds[j][jb].jtag == itag) break; - if (jb == jbonds) - error->one(FLERR,"Fix hyper/local could not find duplicate bond"); - if (bonds[i][m].boostcoeff != bonds[j][jb].boostcoeff) - checkcoeff_count++; - } - } - } } /* ---------------------------------------------------------------------- */ @@ -837,159 +782,215 @@ void FixHyperLocal::build_bond_list(int natom) nevent_atom += natom; } - // trigger Dcut neighbor list build - // NOTE: turn off special bonds in this Dcut neigh list? + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; - neighbor->build_one(list); + // acquire old bond coeffs so can persist them in new blist + // while loop is to allow new value of maxcoeffperatom + // will loop at most 2 times, just once when maxcoeffperatom is large enough + // just reverse comm needed, + // b/c new bond list will be bonds of current owned atoms - // make copy of old bonds to preserve boostcoeffs for bonds that persist - // allocate new numbond + tagint *tag = atom->tag; - OneBond **old_bonds = bonds; - int *old_numbond = numbond; + if (maxcoeff < nall) { + maxcoeff = atom->nmax; + grow_coeff(); + } - int nmax = atom->nmax; - memory->create(numbond,nmax,"hyper/local:numbond"); + while (1) { + if (firstflag) break; + for (i = 0; i < nall; i++) numcoeff[i] = 0; - // old_nall = value of nall at time bonds are built - // reallocate new xold and tagold if necessary + for (m = 0; m < nblocal; m++) { + i = blist[m].i; + j = blist[m].j; + + if (numcoeff[i] < maxcoeffperatom) { + clist[i][numcoeff[i]].biascoeff = blist[m].biascoeff; + clist[i][numcoeff[i]].jtag = tag[j]; + } + numcoeff[i]++; + + if (numcoeff[j] < maxcoeffperatom) { + clist[j][numcoeff[j]].biascoeff = blist[m].biascoeff; + clist[j][numcoeff[i]].jtag = tag[i]; + } + numcoeff[j]++; + } + + int maxcol = 0; + for (i = 0; i < nall; i++) maxcol = MAX(maxcol,numcoeff[i]); + int maxcolall; + MPI_Allreduce(&maxcol,&maxcolall,1,MPI_INT,MPI_MAX,world); + + if (maxcolall > maxcoeffperatom) { + maxcoeffperatom = maxcolall; + grow_coeff(); + memory->destroy(clist); + memory->create(clist,maxcoeff,maxcoeffperatom,"hyper/local:clist"); + continue; + } + + commflag = BIASCOEFF; + comm->reverse_comm_fix(this); + + maxcol = 0; + for (i = 0; i < nall; i++) maxcol = MAX(maxcol,numcoeff[i]); + MPI_Allreduce(&maxcol,&maxcolall,1,MPI_INT,MPI_MAX,world); + if (maxcolall <= maxcoeffperatom) break; + + maxcoeffperatom = maxcolall; + grow_coeff(); + } + + // reallocate vectors that are maxnew xold and tagold if necessary // initialize xold to current coords // initialize tagold to zero, so atoms not in neighbor list will remain zero - old_nall = atom->nlocal + atom->nghost; - - if (old_nall > maxold) { + if (nlocal > maxlocal) { + memory->destroy(eligible); + memory->destroy(numbond); + memory->destroy(maxhalf); + maxlocal = nlocal; + memory->create(eligible,maxlocal,"hyper/local:eligible"); + memory->create(numbond,maxlocal,"hyper/local:numbond"); + memory->create(maxhalf,maxlocal,"hyper/local:maxhalf"); + } + + if (nall > maxall) { memory->destroy(xold); memory->destroy(tagold); memory->destroy(old2now); - maxold = atom->nmax; - memory->create(xold,maxold,3,"hyper/local:xold"); - memory->create(tagold,maxold,"hyper/local:tagold"); - memory->create(old2now,maxold,"hyper/local:old2now"); + maxall = atom->nmax; + memory->create(xold,maxall,3,"hyper/local:xold"); + memory->create(tagold,maxall,"hyper/local:tagold"); + memory->create(old2now,maxall,"hyper/local:old2now"); } + // nlocal_old = value of nlocal at time bonds are built + // nall_old = value of nall at time bonds are built + // archive current peratom info in old vecs + + nlocal_old = nlocal; + nall_old = nall; + double **x = atom->x; - memcpy(&xold[0][0],&x[0][0],3*old_nall*sizeof(double)); - for (i = 0; i < old_nall; i++) tagold[i] = 0; + memcpy(&xold[0][0],&x[0][0],3*nall*sizeof(double)); + for (i = 0; i < nall; i++) tagold[i] = 0; + for (i = 0; i < nlocal; i++) numbond[i] = 0; - // create and populate new bonds data struct - // while loop allows maxbondperatom to increase once if necessary - // don't know new maxbondperatom value until end of loop - // in practice maxbondperatom will hardly ever increase - // since there is a physical max value + // trigger builds for both neighbor lists + // NOTE: insure the I atoms are in same order? + + neighbor->build_one(listfull); + neighbor->build_one(listhalf); - tagint *tag = atom->tag; - int *mask = atom->mask; - int nlocal = atom->nlocal; + // set tagold = 1 for all J atoms used in full neighbor list + // tagold remains 0 for unused atoms, skipped in pre_neighbor - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; + inum = listfull->inum; + ilist = listfull->ilist; + numneigh = listfull->numneigh; + firstneigh = listfull->firstneigh; - while (1) { - bonds = (OneBond **) memory->create(bonds,nmax,maxbondperatom, - "hyper/local:bonds"); - if (bonds) memset(bonds[0],0,nmax*sizeof(OneBond)); - for (i = 0; i < nlocal; i++) numbond[i] = 0; - - // identify bonds assigned to each owned atom - // do not create a bond between two non-group atoms - // set tagold = global ID for all I,J atoms used in neighbor list - // tagold remains 0 for unused atoms, skipped in pre_neighbor - - int nbondmax = 0; - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itag = tag[i]; - tagold[i] = tag[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - nbond = 0; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - jtag = tag[j]; - tagold[j] = jtag; - - // skip if neither atom I or J are in fix group - // order IJ to insure IJ and JI bonds are stored consistently - - if (!(mask[i] & groupbit) && !(mask[j] & groupbit)) continue; - - if (itag < jtag) { - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - } else { - delx = x[j][0] - xtmp; - dely = x[j][1] - ytmp; - delz = x[j][2] - ztmp; - } - - rsq = delx*delx + dely*dely + delz*delz; - - // NOTE: could create two bonds for IJ both owned from one calc? - // have to skip one of 2 bonds in that case - - if (rsq < cutbondsq) { - if (nbond >= maxbondperatom) { - nbond++; - continue; - } - - bonds[i][nbond].r0 = sqrt(rsq); - bonds[i][nbond].jtag = tag[j]; - bonds[i][nbond].j = j; - - if (firstflag) oldcoeff = 0.0; - else { - oldcoeff = 0.0; - jtag = tag[j]; - n = old_numbond[i]; - for (m = 0; m < n; m++) { - if (old_bonds[i][m].jtag == jtag) { - oldcoeff = old_bonds[i][m].boostcoeff; - break; - } - } - } - - if (oldcoeff > 0.0) bonds[i][nbond].boostcoeff = oldcoeff; - else { - bonds[i][nbond].boostcoeff = BOOSTINIT; - nnewbond++; - } - nbond++; - } - } - numbond[i] = nbond; - nbondmax = MAX(nbondmax,nbond); + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + tagold[j] = 1; } - - // maxbondperatom must increase uniformly on all procs - // since bonds are comunicated when atoms migrate - - int allnbondmax; - MPI_Allreduce(&nbondmax,&allnbondmax,1,MPI_INT,MPI_MAX,world); - if (allnbondmax <= maxbondperatom) break; - - maxbondperatom = allnbondmax; - memory->destroy(bonds); } - // deallocate old_bonds and old_numbond + // identify bonds assigned to each owned atom + // do not create a bond between two non-group atoms - memory->destroy(old_bonds); - memory->destroy(old_numbond); + int *mask = atom->mask; + + inum = listhalf->inum; + ilist = listhalf->ilist; + numneigh = listhalf->numneigh; + firstneigh = listhalf->firstneigh; + + bigint bondcount = 0; + nblocal = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + itag = tag[i]; + tagold[i] = tag[i]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + nbond = 0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + jtag = tag[j]; + tagold[j] = jtag; + + // skip if neither atom I or J are in fix group + + if (!(mask[i] & groupbit) && !(mask[j] & groupbit)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cutbondsq) { + nbond++; + + if (nblocal == maxbond) grow_bond(); + blist[nblocal].i = i; + blist[nblocal].j = j; + blist[nblocal].iold = i; + blist[nblocal].jold = j; + blist[nblocal].r0 = sqrt(rsq); + + // set biascoeff to old coeff for same I,J pair or to default + + if (firstflag) oldcoeff = 0.0; + else { + oldcoeff = 0.0; + jtag = tag[j]; + n = numcoeff[i]; + for (m = 0; m < n; m++) { + if (clist[i][m].jtag == jtag) { + oldcoeff = clist[i][m].biascoeff; + break; + } + } + } + + if (oldcoeff > 0.0) blist[nblocal].biascoeff = oldcoeff; + else { + blist[nblocal].biascoeff = COEFFINIT; + nnewbond++; + } + + nblocal++; + } + } + + numbond[i] = nbond; + bondcount += nbond; + } + + MPI_Allreduce(&bondcount,&allbonds,1,MPI_LMP_BIGINT,MPI_SUM,world); time2 = MPI_Wtime(); + if (firstflag) nnewbond = 0; else { time_bondbuild += time2-time1; @@ -1016,13 +1017,13 @@ int FixHyperLocal::pack_forward_comm(int n, int *list, double *buf, buf[m++] = maxstrain[j]; } - // STRAINREGION - // pack maxstrain_region vector + // STRAINDOMAIN + // pack maxstrain_domain vector - } else if (commflag == STRAINREGION) { + } else if (commflag == STRAINDOMAIN) { for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = maxstrain_region[j]; + buf[m++] = maxstrain_domain[j]; } // BIASFLAG @@ -1056,11 +1057,11 @@ void FixHyperLocal::unpack_forward_comm(int n, int first, double *buf) } // STRAINREGION - // unpack maxstrain_region vector + // unpack maxstrain_domain vector - } else if (commflag == STRAINREGION) { + } else if (commflag == STRAINDOMAIN) { for (i = first; i < last; i++) { - maxstrain_region[i] = buf[m++]; + maxstrain_domain[i] = buf[m++]; } // BIASFLAG @@ -1077,17 +1078,25 @@ void FixHyperLocal::unpack_forward_comm(int n, int first, double *buf) int FixHyperLocal::pack_reverse_comm(int n, int first, double *buf) { - int i,m,last; + int i,j,m,last; m = 0; last = first + n; - // STRAINREGION - // pack maxstrain_region vector + // STRAIN + // pack maxstrain vector - if (commflag == STRAINREGION) { + if (commflag == STRAIN) { for (i = first; i < last; i++) { - buf[m++] = maxstrain_region[i]; + buf[m++] = maxstrain[i]; + } + + // STRAINDOMAIN + // pack maxstrain_domain vector + + } else if (commflag == STRAINDOMAIN) { + for (i = first; i < last; i++) { + buf[m++] = maxstrain_domain[i]; } // BIASFLAG @@ -1097,6 +1106,20 @@ int FixHyperLocal::pack_reverse_comm(int n, int first, double *buf) for (i = first; i < last; i++) { buf[m++] = ubuf(biasflag[i]).d; } + + // BIASCOEFF + // pack list of biascoeffs + + } else if (commflag == BIASCOEFF) { + int ncoeff; + for (i = first; i < last; i++) { + ncoeff = numcoeff[i]; + buf[m++] = ubuf(ncoeff).d; + for (j = 0; j < ncoeff; j++) { + buf[m++] = clist[i][j].biascoeff; + buf[m++] = ubuf(clist[i][j].jtag).d; + } + } } return m; @@ -1106,17 +1129,29 @@ int FixHyperLocal::pack_reverse_comm(int n, int first, double *buf) void FixHyperLocal::unpack_reverse_comm(int n, int *list, double *buf) { - int i,j,m; + int i,j,k,m; m = 0; - // STRAINREGION - // unpack maxstrain_region vector - - if (commflag == STRAINREGION) { + // STRAIN + // unpack maxstrain vector + // use MAX, b/c want maximum abs value strain for each atom's bonds + + if (commflag == STRAIN) { for (i = 0; i < n; i++) { j = list[i]; - maxstrain_region[j] += buf[m++]; + maxstrain[j] = MAX(maxstrain[j],buf[m]); + m++; + } + + // STRAINDOMAIN + // unpack maxstrain_domain vector + // use SUM, b/c exactly one ghost or owned value per atom ID is non-zero + + } else if (commflag == STRAINDOMAIN) { + for (i = 0; i < n; i++) { + j = list[i]; + maxstrain_domain[j] += buf[m++]; } // BIASFLAG @@ -1127,84 +1162,49 @@ void FixHyperLocal::unpack_reverse_comm(int n, int *list, double *buf) j = list[i]; biasflag[j] = (tagint) ubuf(buf[m++]).i; } + + // BIASCOEFF + // unpack list of biascoeffs and add to atom J's list + // protect against overflow of clist columns + // if that happens, caller will realloc clist and reverse comm again + + } else if (commflag == BIASFLAG) { + int ncoeff; + for (i = 0; i < n; i++) { + j = list[i]; + ncoeff = (int) ubuf(buf[m++]).i; + for (k = 0; k < ncoeff; k++) { + if (numcoeff[j] < maxcoeffperatom) { + clist[j][numcoeff[j]].biascoeff = buf[m++]; + clist[j][numcoeff[j]].jtag = (tagint) ubuf(buf[m++]).i; + } else m += 2; + numcoeff[j]++; + } + } } } /* ---------------------------------------------------------------------- - allocate atom-based arrays + grow bond list by a chunk ------------------------------------------------------------------------- */ -void FixHyperLocal::grow_arrays(int nmax) +void FixHyperLocal::grow_bond() { - // NOTE: not all of these need to be Nmax in length, could allocate elsewhere - - memory->grow(maxstrain,nmax,"hyper/local:maxstrain"); - memory->grow(maxstrain_bondindex,nmax,"hyper/local:maxstrain_bondindex"); - memory->grow(maxstrain_region,nmax,"hyper/local:maxstrain_region"); - if (checkbias) memory->grow(biasflag,nmax,"hyper/local:biasflag"); - - memory->grow(numbond,nmax,"hyper/local:numbond"); - memory->grow(bonds,nmax,maxbondperatom,"hyper/local:bonds"); - - // zero so valgrind does not complain about memcpy() in copy() - // also so loops in pre_neighbor() are OK before - // bonds are setup for the first time - - if (bonds) { - memset(bonds[maxbond],0,(nmax-maxbond)*maxbondperatom*sizeof(OneBond)); - memset(&numbond[maxbond],0,(nmax-maxbond)*sizeof(int)); - maxbond = nmax; - } + if (maxbond + DELTABOND > MAXSMALLINT) + error->one(FLERR,"Fix hyper/local bond count is too big"); + maxbond += DELTABOND; + blist = (OneBond *) + memory->srealloc(blist,maxbond*sizeof(OneBond),"hyper/local:blist"); } /* ---------------------------------------------------------------------- - copy values within local atom-based arrays + reallocate 2-dimensional clist ------------------------------------------------------------------------- */ -void FixHyperLocal::copy_arrays(int i, int j, int /* delflag */) +void FixHyperLocal::grow_coeff() { - // avoid valgrind copy-to-self warning - - if (i != j) memcpy(bonds[j],bonds[i],numbond[i]*sizeof(OneBond)); - numbond[j] = numbond[i]; -} - -/* ---------------------------------------------------------------------- - pack values in local atom-based array for exchange with another proc -------------------------------------------------------------------------- */ - -int FixHyperLocal::pack_exchange(int i, double *buf) -{ - int m = 1; - int n = numbond[i]; - buf[m++] = ubuf(n).d; - for (int j = 0; j < n; j++) { - buf[m++] = bonds[i][j].r0; - buf[m++] = bonds[i][j].boostcoeff; - buf[m++] = ubuf(bonds[i][j].jtag).d; - buf[m++] = ubuf(bonds[i][j].j).d; - } - - buf[0] = m; - return m; -} - -/* ---------------------------------------------------------------------- - unpack values in local atom-based array from exchange with another proc -------------------------------------------------------------------------- */ - -int FixHyperLocal::unpack_exchange(int nlocal, double *buf) -{ - int m = 1; - int n = numbond[nlocal] = (int) ubuf(buf[m++]).i; - for (int j = 0; j < n; j++) { - bonds[nlocal][j].r0 = buf[m++]; - bonds[nlocal][j].boostcoeff = buf[m++]; - bonds[nlocal][j].jtag = (tagint) ubuf(buf[m++]).i; - bonds[nlocal][j].j = (int) ubuf(buf[m++]).i; - } - - return m; + memory->destroy(clist); + memory->create(clist,maxcoeff,maxcoeffperatom,"hyper/local:clist"); } /* ---------------------------------------------------------------------- */ @@ -1222,7 +1222,7 @@ double FixHyperLocal::compute_vector(int i) { // 23 vector outputs returned for i = 0-22 - // i = 0 = # of boosted bonds on this step + // i = 0 = # of biased bonds on this step // i = 1 = max strain of any bond on this step // i = 2 = average bias potential for all bonds on this step // i = 3 = ave bonds/atom on this step @@ -1231,7 +1231,7 @@ double FixHyperLocal::compute_vector(int i) // i = 5 = fraction of steps and bonds with no bias during this run // i = 6 = max drift distance of any atom during this run // i = 7 = max bond length during this run - // i = 8 = average # of boosted bonds/step during this run + // i = 8 = average # of biased bonds/step during this run // i = 9 = average bias potential for all bonds during this run // i = 10 = max bias potential for any bond during this run // i = 11 = min bias potential for any bond during this run @@ -1241,9 +1241,14 @@ double FixHyperLocal::compute_vector(int i) // any maxstrain during this run // i = 14 = count of ghost atoms that could not be found // by any proc at any reneighbor step during this run + + // NOTE: these 2 are no longer relevant // i = 15 = count of lost bond partners during this run // i = 16 = average bias coeff for lost bond partners during this run + // i = 17 = count of bias overlaps found during this run + + // NOTE: this is no longer relevant // i = 18 = count of non-matching bias coefficients found during this run // i = 19 = cummulative hyper time @@ -1252,9 +1257,9 @@ double FixHyperLocal::compute_vector(int i) // i = 22 = cummulative # of new bonds formed since fix created if (i == 0) { - int nboostall; - MPI_Allreduce(&nboost,&nboostall,1,MPI_INT,MPI_SUM,world); - return (double) nboostall; + int nbiasall; + MPI_Allreduce(&nbias,&nbiasall,1,MPI_INT,MPI_SUM,world); + return (double) nbiasall; } if (i == 1) { @@ -1269,33 +1274,30 @@ double FixHyperLocal::compute_vector(int i) } if (i == 2) { - if (allboost && allbonds) return allboost/allbonds * vmax; + if (allbias && allbonds) return allbias/allbonds * vmax; return 1.0; } - if (i == 3) return 1.0*allbonds/atom->natoms; + if (i == 3) { + return 2.0*allbonds/atom->natoms; + } if (i == 4) { const int nlocal = atom->nlocal; - bigint nbonds = 0; - for (int j = 0; j < nlocal; j++) - nbonds += numbond[j]; - bigint allbonds; - MPI_Allreduce(&nbonds,&allbonds,1,MPI_LMP_BIGINT,MPI_SUM,world); bigint allneigh,thisneigh; - thisneigh = list->ipage->ndatum; + thisneigh = listfull->ipage->ndatum; MPI_Allreduce(&thisneigh,&allneigh,1,MPI_LMP_BIGINT,MPI_SUM,world); const double natoms = atom->natoms; const double neighsperatom = static_cast(allneigh)/natoms; - const double bondsperatom = 0.5*static_cast(allbonds)/natoms; + const double bondsperatom = static_cast(allbonds)/natoms; return neighsperatom * bondsperatom; } if (i == 5) { - int allboost_running,allnobias_running; - MPI_Allreduce(&nboost_running,&allboost_running,1,MPI_INT,MPI_SUM,world); + int allbias_running,allnobias_running; + MPI_Allreduce(&nbias_running,&allbias_running,1,MPI_INT,MPI_SUM,world); MPI_Allreduce(&nobias_running,&allnobias_running,1,MPI_INT,MPI_SUM,world); - if (allboost_running) return 1.0*allnobias_running / allboost_running; + if (allbias_running) return 1.0*allnobias_running / allbias_running; return 0.0; } @@ -1313,26 +1315,26 @@ double FixHyperLocal::compute_vector(int i) if (i == 8) { if (update->ntimestep == update->firststep) return 0.0; - int allboost_running; - MPI_Allreduce(&nboost_running,&allboost_running,1,MPI_INT,MPI_SUM,world); - return 1.0*allboost_running / (update->ntimestep - update->firststep); + int allbias_running; + MPI_Allreduce(&nbias_running,&allbias_running,1,MPI_INT,MPI_SUM,world); + return 1.0*allbias_running / (update->ntimestep - update->firststep); } if (i == 9) { if (update->ntimestep == update->firststep) return 0.0; - return sumboostcoeff * vmax / (update->ntimestep - update->firststep); + return sumbiascoeff * vmax / (update->ntimestep - update->firststep); } if (i == 10) { - double allboostcoeff; - MPI_Allreduce(&maxboostcoeff,&allboostcoeff,1,MPI_DOUBLE,MPI_MAX,world); - return allboostcoeff * vmax; + double allbiascoeff; + MPI_Allreduce(&maxbiascoeff,&allbiascoeff,1,MPI_DOUBLE,MPI_MAX,world); + return allbiascoeff * vmax; } if (i == 11) { - double allboostcoeff; - MPI_Allreduce(&minboostcoeff,&allboostcoeff,1,MPI_DOUBLE,MPI_MAX,world); - return allboostcoeff * vmax; + double allbiascoeff; + MPI_Allreduce(&minbiascoeff,&allbiascoeff,1,MPI_DOUBLE,MPI_MAX,world); + return allbiascoeff * vmax; } if (i == 12) return rmaxever; @@ -1344,35 +1346,14 @@ double FixHyperLocal::compute_vector(int i) return 1.0*allghost_toofar; } - if (i == 15) { - int alllost; - MPI_Allreduce(&lostbond_partner,&alllost,1,MPI_INT,MPI_SUM,world); - return 1.0*alllost; - } - - if (i == 16) { - int alllost; - MPI_Allreduce(&lostbond_partner,&alllost,1,MPI_INT,MPI_SUM,world); - double allcoeff; - MPI_Allreduce(&lostbond_coeff,&allcoeff,1,MPI_DOUBLE,MPI_SUM,world); - if (alllost == 0) return 0; - return allcoeff/alllost; - } - if (i == 17) { int allclose; MPI_Allreduce(&checkbias_count,&allclose,1,MPI_INT,MPI_SUM,world); return 1.0*allclose; } - if (i == 18) { - int allcoeff; - MPI_Allreduce(&checkcoeff_count,&allcoeff,1,MPI_INT,MPI_SUM,world); - return 1.0*allcoeff; - } - if (i == 19) { - return boosttarget * update->dt * (update->ntimestep - starttime); + return boost_target * update->dt * (update->ntimestep - starttime); } if (i == 20) return (double) nevent; @@ -1406,20 +1387,20 @@ double FixHyperLocal::query(int i) if (i == 8) return compute_vector(22); // number of new bonds if (i == 9) return 1.0*maxbondperatom; // max bonds/atom - if (i == 10) return compute_vector(8); // ave # of boosted bonds/step - if (i == 11) return compute_vector(9); // ave boost coeff over all bonds - if (i == 12) return compute_vector(10); // max boost cooef for any bond - if (i == 13) return compute_vector(11); // max boost cooef for any bond + if (i == 10) return compute_vector(8); // ave # of biased bonds/step + if (i == 11) return compute_vector(9); // ave bias coeff over all bonds + if (i == 12) return compute_vector(10); // max bias cooef for any bond + if (i == 13) return compute_vector(11); // max bias cooef for any bond if (i == 14) return compute_vector(4); // neighbor bonds/bond - if (i == 15) return compute_vector(2); // ave boost cooef now + if (i == 15) return compute_vector(2); // ave bias coeff now if (i == 16) return time_bondbuild; // CPU time for bond_build calls if (i == 17) return rmaxever; // ghost atom distance for < maxstrain if (i == 18) return rmaxeverbig; // ghost atom distance for any strain if (i == 19) return compute_vector(14); // count of ghost atoms not found - if (i == 20) return compute_vector(15); // count of lost bond partners - if (i == 21) return compute_vector(16); // ave bias coeff of long bonds + //if (i == 20) return compute_vector(15); // count of lost bond partners + //if (i == 21) return compute_vector(16); // ave bias coeff of long bonds if (i == 22) return compute_vector(17); // count of bias overlaps - if (i == 23) return compute_vector(18); // count of non-matching bias coeffs + //if (i == 23) return compute_vector(18); // count of non-matching bias coeffs error->all(FLERR,"Invalid query to fix hyper/local"); @@ -1433,11 +1414,15 @@ double FixHyperLocal::query(int i) double FixHyperLocal::memory_usage() { int nmax = atom->nmax; - double bytes = 2*nmax * sizeof(int); // numbond, maxstrain_bondindex - bytes = 2*nmax * sizeof(double); // maxstrain, maxstrain_region - bytes += maxbondperatom*nmax * sizeof(OneBond); // bonds - bytes += 3*maxold * sizeof(double); // xold - bytes += maxold * sizeof(tagint); // tagold - bytes += maxold * sizeof(int); // old2now + double bytes = maxbond * sizeof(OneBond); // bond list + bytes += 3*maxlocal * sizeof(int); // numbond,maxhalf,eligible + bytes += maxlocal * sizeof(double); // maxhalfstrain + bytes += maxall * sizeof(int); // old2now + bytes += maxall * sizeof(tagint); // tagold + bytes += 3*maxall * sizeof(double); // xold + bytes += 2*nmax * sizeof(double); // maxstrain,maxstrain_domain + if (checkbias) bytes += nmax * sizeof(tagint); // biasflag + bytes += maxcoeff*maxcoeffperatom * sizeof(OneCoeff); // clist + bytes += maxcoeff * sizeof(int); // numcoeff return bytes; } diff --git a/src/REPLICA/fix_hyper_local.h b/src/REPLICA/fix_hyper_local.h index 147e3ef1ef..67361ce3ac 100644 --- a/src/REPLICA/fix_hyper_local.h +++ b/src/REPLICA/fix_hyper_local.h @@ -45,11 +45,6 @@ class FixHyperLocal : public FixHyper { 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 @@ -62,24 +57,20 @@ class FixHyperLocal : public FixHyper { 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 @@ -91,65 +82,88 @@ class FixHyperLocal : public FixHyper { 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; + 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 + + // list of my owned bonds + // persists on a proc from one event until the next + + 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 + }; + + struct OneBond *blist; // list of owned bonds + int nblocal; // # of owned bonds + int maxbond; // allocated size of blist + + // 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 + + 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 *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 + + 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 + + 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_domain; // max-strain of any neighbor atom J of atom I + // for owned and ghost atoms + tagint *biasflag; // atoms in biased bonds marked with bond partner + // for owned and ghost atoms + + // data struct used to persist biascoeffs when bond list is re-created + + 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; - // 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 + // private methods - class NeighList *list; // full neigh list up to Dcut distance - // created only when bonds are rebuilt - - 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 - - double **xold; // coords of owned+ghost atoms when bonds created - tagint *tagold; // global IDs of owned+ghost atoms when b created - - int maxold; // allocated size of old2now - int maxbond; // allocated size of bonds - int old_nall; // nlocal+nghost when old2now was last setup - - 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 - }; - - struct OneBond **bonds; // 2d array of bonds for owned atoms - int *numbond; // number of bonds for each owned atom - - 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 - // 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 - - int maxboost; // allocated size of boost list - int nboost; // # of boosted bonds I own - int *boost; // index of atom I in each boosted bond - - // histogramming of bond boost cooeficients - - 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(); }; }