diff --git a/src/USER-ACKLAND/compute_ackland_atom.cpp b/src/USER-ACKLAND/compute_ackland_atom.cpp
index 6ded5cf4f..536631f31 100644
--- a/src/USER-ACKLAND/compute_ackland_atom.cpp
+++ b/src/USER-ACKLAND/compute_ackland_atom.cpp
@@ -1,385 +1,392 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: G. Ziegenhain, gerolf@ziegenhain.com
Copyright (C) 2007
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_ackland_atom.h"
#include "atom.h"
#include "modify.h"
#include "update.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "pair.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include <math.h>
using namespace LAMMPS_NS;
enum{UNKNOWN,BCC,FCC,HCP,ICO};
/* ---------------------------------------------------------------------- */
ComputeAcklandAtom::ComputeAcklandAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all("Illegal compute ackland/atom command");
peratom_flag = 1;
size_peratom = 0;
nmax = 0;
structure = NULL;
maxneigh = 0;
distsq = NULL;
nearest = NULL;
nearest_n0 = NULL;
nearest_n1 = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeAcklandAtom::~ComputeAcklandAtom()
{
memory->sfree(structure);
memory->sfree(distsq);
memory->sfree(nearest);
memory->sfree(nearest_n0);
memory->sfree(nearest_n1);
}
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::init()
{
// need an occasional full neighbor list
int irequest = neighbor->request((void *) this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->occasional = 1;
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"ackland/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning("More than one compute ackland/atom");
}
+/* ---------------------------------------------------------------------- */
+
+void ComputeAcklandAtom::init_list(int id, NeighList *ptr)
+{
+ list = ptr;
+}
+
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::compute_peratom()
{
int i,j,ii,jj,k,n,inum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq,value;
int *ilist,*jlist,*numneigh,**firstneigh;
double pairs[66];
int chi[8];
// grow structure array if necessary
if (atom->nlocal > nmax) {
memory->sfree(structure);
nmax = atom->nmax;
structure = (double *)
memory->smalloc(nmax*sizeof(double),"compute/ackland/atom:ackland");
scalar_atom = structure;
}
// invoke half neighbor list (will copy or build if necessary)
neighbor->build_one(list->index);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// compute structure parameter for each atom in group
// use full neighbor list
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
double cutsq = force->pair->cutforce * force->pair->cutforce;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
// ensure distsq and nearest arrays are long enough
if (jnum > maxneigh) {
memory->sfree(distsq);
memory->sfree(nearest);
memory->sfree(nearest_n0);
memory->sfree(nearest_n1);
maxneigh = jnum;
distsq = (double *) memory->smalloc(maxneigh*sizeof(double),
"compute/ackland/atom:distsq");
nearest = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/ackland/atom:nearest");
nearest_n0 = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/ackland/atom:nearest_n0");
nearest_n1 = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/ackland/atom:nearest_n1");
}
// loop over list of all neighbors within force cutoff
// distsq[] = distance sq to each
// nearest[] = atom indices of neighbors
n = 0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (j >= nall) j %= nall;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutsq) {
distsq[n] = rsq;
nearest[n++] = j;
}
}
// Select 6 nearest neighbors
select2(6,n,distsq,nearest);
// Mean squared separation
double r0_sq = 0.;
for (j = 0; j < 6; j++)
r0_sq += distsq[j];
r0_sq /= 6.;
// n0 near neighbors with: distsq<1.45*r0_sq
// n1 near neighbors with: distsq<1.55*r0_sq
double n0_dist_sq = 1.45*r0_sq,
n1_dist_sq = 1.55*r0_sq;
int n0 = 0, n1 = 0;
for (j = 0; j < n; j++) {
if (distsq[j] < n1_dist_sq) {
nearest_n1[n1++] = nearest[j];
if (distsq[j] < n0_dist_sq) {
nearest_n0[n0++] = nearest[j];
}
}
}
// Evaluate all angles <(r_ij,rik) forall n0 particles with: distsq<1.45*r0_sq
double bond_angle;
double norm_j, norm_k;
chi[0] = chi[1] = chi[2] = chi[3] = chi[4] = chi[5] = chi[6] = chi[7] = 0;
double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik;
for (j = 0; j < n0; j++) {
x_ij = x[i][0]-x[nearest_n0[j]][0];
y_ij = x[i][1]-x[nearest_n0[j]][1];
z_ij = x[i][2]-x[nearest_n0[j]][2];
norm_j = sqrt (x_ij*x_ij + y_ij*y_ij + z_ij*z_ij);
if (norm_j <= 0.) continue;
for (k = j+1; k < n0; k++) {
x_ik = x[i][0]-x[nearest_n0[k]][0];
y_ik = x[i][1]-x[nearest_n0[k]][1];
z_ik = x[i][2]-x[nearest_n0[k]][2];
norm_k = sqrt (x_ik*x_ik + y_ik*y_ik + z_ik*z_ik);
if (norm_k <= 0.)
continue;
bond_angle = (x_ij*x_ik + y_ij*y_ik + z_ij*z_ik) / (norm_j*norm_k);
// Histogram for identifying the relevant peaks
if (-1. <= bond_angle && bond_angle < -0.945) { chi[0]++; }
else if (-0.945 <= bond_angle && bond_angle < -0.915) { chi[1]++; }
else if (-0.915 <= bond_angle && bond_angle < -0.755) { chi[2]++; }
else if (-0.755 <= bond_angle && bond_angle < -0.195) { chi[3]++; }
else if (-0.195 <= bond_angle && bond_angle < 0.195) { chi[4]++; }
else if (0.195 <= bond_angle && bond_angle < 0.245) { chi[5]++; }
else if (0.245 <= bond_angle && bond_angle < 0.795) { chi[6]++; }
else if (0.795 <= bond_angle && bond_angle < 1.) { chi[7]++; }
}
}
// Deviations from the different lattice structures
double delta_bcc = 0.35*chi[4]/(double)(chi[5]+chi[6]-chi[4]),
delta_cp = fabs(1.-(double)chi[6]/24.),
delta_fcc = 0.61*(fabs((double)(chi[0]+chi[1]-6.))+(double)chi[2])/6.,
delta_hcp = (fabs((double)chi[0]-3.)+fabs((double)chi[0]+(double)chi[1]+(double)chi[2]+(double)chi[3]-9.))/12.;
// Identification of the local structure according to the reference
if (chi[0] == 7) { delta_bcc = 0.; }
else if (chi[0] == 6) { delta_fcc = 0.; }
else if (chi[0] <= 3) { delta_hcp = 0.; }
if (chi[7] > 0.)
structure[i] = UNKNOWN;
else
if (chi[4] < 3.)
{
if (n1 > 13 || n1 < 11)
structure[i] = UNKNOWN;
else
structure[i] = ICO;
} else
if (delta_bcc <= delta_cp)
{
if (n1 < 11)
structure[i] = UNKNOWN;
else
structure[i] = BCC;
} else
if (n1 > 12 || n1 < 11)
structure[i] = UNKNOWN;
else
if (delta_fcc < delta_hcp)
structure[i] = FCC;
else
structure[i] = HCP;
} // end loop over all particles
}
}
/* ----------------------------------------------------------------------
2 select routines from Numerical Recipes (slightly modified)
find k smallest values in array of length n
2nd routine sorts auxiliary array at same time
------------------------------------------------------------------------- */
#define SWAP(a,b) tmp = a; a = b; b = tmp;
#define ISWAP(a,b) itmp = a; a = b; b = itmp;
void ComputeAcklandAtom::select(int k, int n, double *arr)
{
int i,ir,j,l,mid;
double a,tmp;
arr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir])
}
return;
} else {
mid=(l+ir) >> 1;
SWAP(arr[mid],arr[l+1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
}
i = l+1;
j = ir;
a = arr[l+1];
for (;;) {
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
if (j < i) break;
SWAP(arr[i],arr[j])
}
arr[l+1] = arr[j];
arr[j] = a;
if (j >= k) ir = j-1;
if (j <= k) l = i;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::select2(int k, int n, double *arr, int *iarr)
{
int i,ir,j,l,mid,ia,itmp;
double a,tmp;
arr--;
iarr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir])
ISWAP(iarr[l],iarr[ir])
}
return;
} else {
mid=(l+ir) >> 1;
SWAP(arr[mid],arr[l+1])
ISWAP(iarr[mid],iarr[l+1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
ISWAP(iarr[l],iarr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
ISWAP(iarr[l+1],iarr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
ISWAP(iarr[l],iarr[l+1])
}
i = l+1;
j = ir;
a = arr[l+1];
ia = iarr[l+1];
for (;;) {
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
if (j < i) break;
SWAP(arr[i],arr[j])
ISWAP(iarr[i],iarr[j])
}
arr[l+1] = arr[j];
arr[j] = a;
iarr[l+1] = iarr[j];
iarr[j] = ia;
if (j >= k) ir = j-1;
if (j <= k) l = i;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeAcklandAtom::memory_usage()
{
double bytes = nmax * sizeof(double);
return bytes;
}
diff --git a/src/USER-ACKLAND/compute_ackland_atom.h b/src/USER-ACKLAND/compute_ackland_atom.h
index a4d56713f..1c10110a4 100644
--- a/src/USER-ACKLAND/compute_ackland_atom.h
+++ b/src/USER-ACKLAND/compute_ackland_atom.h
@@ -1,42 +1,43 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef COMPUTE_ACKLAND_ATOM_H
#define COMPUTE_ACKLAND_ATOM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeAcklandAtom : public Compute {
public:
ComputeAcklandAtom(class LAMMPS *, int, char **);
~ComputeAcklandAtom();
void init();
+ void init_list(int, class NeighList *);
void compute_peratom();
double memory_usage();
private:
int nmax,maxneigh;
double *distsq;
int *nearest, *nearest_n0, *nearest_n1;
double *structure;
class NeighList *list;
void select(int, int, double *);
void select2(int, int, double *, int *);
};
}
#endif
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index fa4790ba2..6b906ed99 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -1,1479 +1,1479 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author (triclinic and multi-neigh) : Pieter in 't Veld (SNL)
------------------------------------------------------------------------- */
#include "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "limits.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "group.h"
#include "modify.h"
#include "fix.h"
#include "compute.h"
#include "update.h"
#include "respa.h"
#include "output.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define RQDELTA 1
#define EXDELTA 1
#define LB_FACTOR 1.5
#define SMALL 1.0e-6
#define BIG 1.0e20
#define CUT2BIN_RATIO 100
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
enum{NSQ,BIN,MULTI}; // also in neigh_list.cpp
/* ---------------------------------------------------------------------- */
Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp)
{
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
style = BIN;
every = 1;
delay = 10;
dist_check = 1;
pgsize = 10000;
oneatom = 2000;
cutneighsq = NULL;
cuttype = NULL;
cuttypesq = NULL;
fixchecklist = NULL;
// coords at last neighboring
maxhold = 0;
xhold = NULL;
// binning
maxhead = 0;
binhead = NULL;
maxbin = 0;
bins = NULL;
// pair exclusion list info
nex_type = maxex_type = 0;
ex1_type = ex2_type = NULL;
ex_type = NULL;
nex_group = maxex_group = 0;
ex1_group = ex2_group = ex1_bit = ex2_bit = NULL;
nex_mol = maxex_mol = 0;
ex_mol_group = ex_mol_bit = NULL;
// pair lists
maxlocal = 0;
nlist = 0;
lists = NULL;
pair_build = NULL;
stencil_create = NULL;
blist = glist = slist = NULL;
nrequest = maxrequest = 0;
requests = NULL;
old_style = BIN;
old_nrequest = 0;
old_requests = NULL;
// bond lists
maxbond = 0;
bondlist = NULL;
maxangle = 0;
anglelist = NULL;
maxdihedral = 0;
dihedrallist = NULL;
maximproper = 0;
improperlist = NULL;
}
/* ---------------------------------------------------------------------- */
Neighbor::~Neighbor()
{
memory->destroy_2d_double_array(cutneighsq);
delete [] cuttype;
delete [] cuttypesq;
delete [] fixchecklist;
memory->destroy_2d_double_array(xhold);
memory->sfree(binhead);
memory->sfree(bins);
memory->sfree(ex1_type);
memory->sfree(ex2_type);
memory->destroy_2d_int_array(ex_type);
memory->sfree(ex1_group);
memory->sfree(ex2_group);
delete [] ex1_bit;
delete [] ex2_bit;
memory->sfree(ex_mol_group);
delete [] ex_mol_bit;
for (int i = 0; i < nlist; i++) delete lists[i];
delete [] lists;
delete [] pair_build;
delete [] stencil_create;
delete [] blist;
delete [] glist;
delete [] slist;
for (int i = 0; i < nrequest; i++) delete requests[i];
memory->sfree(requests);
for (int i = 0; i < old_nrequest; i++) delete old_requests[i];
memory->sfree(old_requests);
memory->destroy_2d_int_array(bondlist);
memory->destroy_2d_int_array(anglelist);
memory->destroy_2d_int_array(dihedrallist);
memory->destroy_2d_int_array(improperlist);
}
/* ---------------------------------------------------------------------- */
void Neighbor::init()
{
int i,j,m,n;
ncalls = ndanger = 0;
dimension = domain->dimension;
triclinic = domain->triclinic;
newton_pair = force->newton_pair;
// error check
if (delay > 0 && (delay % every) != 0)
error->all("Neighbor delay must be 0 or multiple of every setting");
// ------------------------------------------------------------------
// settings
// set neighbor cutoffs (force cutoff + skin)
// trigger determines when atoms migrate and neighbor lists are rebuilt
// cutneigh and cutneighsq determine what pairs go into neighbor list
// set to 0 if cutforce = 0
// cutneighmin/max used for neighbor bin sizes
// cutghost determines comm distance = max of cutneigh & skin
// may need ghosts for bonds even if all cutneigh = 0 (pair = NULL)
triggersq = 0.25*skin*skin;
n = atom->ntypes;
if (cutneighsq == NULL) {
cutneighsq = memory->create_2d_double_array(n+1,n+1,"neigh:cutneighsq");
cuttype = new double[n+1];
cuttypesq = new double[n+1];
}
double cutoff,delta,cut;
cutneighmin = BIG;
cutneighmax = 0.0;
for (i = 1; i <= n; i++) {
cuttype[i] = cuttypesq[i] = 0.0;
for (j = 1; j <= n; j++) {
if (force->pair) cutoff = sqrt(force->pair->cutsq[i][j]);
else cutoff = 0.0;
if (cutoff > 0.0) delta = skin;
else delta = 0.0;
cut = cutoff + delta;
cutneighsq[i][j] = cut*cut;
cuttype[i] = MAX(cuttype[i],cut);
cuttypesq[i] = MAX(cuttypesq[i],cut*cut);
cutneighmin = MIN(cutneighmin,cut);
cutneighmax = MAX(cutneighmax,cut);
}
}
cutghost = MAX(cutneighmax,skin);
cutneighmaxsq = cutneighmax * cutneighmax;
// check other classes that can induce reneighboring in decide()
restart_check = 0;
if (output->restart_every) restart_check = 1;
delete [] fixchecklist;
fixchecklist = NULL;
fixchecklist = new int[modify->nfix];
fix_check = 0;
for (i = 0; i < modify->nfix; i++)
if (modify->fix[i]->force_reneighbor)
fixchecklist[fix_check++] = i;
must_check = 0;
if (restart_check || fix_check) must_check = 1;
// set special_flag for 1-2, 1-3, 1-4 neighbors
// flag[0] is not used, flag[1] = 1-2, flag[2] = 1-3, flag[3] = 1-4
// flag = 0 if both LJ/Coulomb special values are 0.0
// flag = 1 if both LJ/Coulomb special values are 1.0
// flag = 2 otherwise or if KSpace solver is enabled
// pairwise portion of KSpace solver uses all 1-2,1-3,1-4 neighbors
if (force->special_lj[1] == 0.0 && force->special_coul[1] == 0.0)
special_flag[1] = 0;
else if (force->special_lj[1] == 1.0 && force->special_coul[1] == 1.0)
special_flag[1] = 1;
else special_flag[1] = 2;
if (force->special_lj[2] == 0.0 && force->special_coul[2] == 0.0)
special_flag[2] = 0;
else if (force->special_lj[2] == 1.0 && force->special_coul[2] == 1.0)
special_flag[2] = 1;
else special_flag[2] = 2;
if (force->special_lj[3] == 0.0 && force->special_coul[3] == 0.0)
special_flag[3] = 0;
else if (force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0)
special_flag[3] = 1;
else special_flag[3] = 2;
if (force->kspace) special_flag[1] = special_flag[2] = special_flag[3] = 2;
// rRESPA cutoffs
int respa = 0;
if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) {
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
}
if (respa) {
double *cut_respa = ((Respa *) update->integrate)->cutoff;
cut_inner_sq = (cut_respa[1] + skin) * (cut_respa[1] + skin);
cut_middle_sq = (cut_respa[3] + skin) * (cut_respa[3] + skin);
cut_middle_inside_sq = (cut_respa[0] - skin) * (cut_respa[0] - skin);
}
// ------------------------------------------------------------------
// xhold, bins, exclusion lists
// free xhold and bins if not needed for this run
if (dist_check == 0) {
memory->destroy_2d_double_array(xhold);
maxhold = 0;
xhold = NULL;
}
if (style == NSQ) {
memory->sfree(bins);
memory->sfree(binhead);
maxbin = maxhead = 0;
binhead = NULL;
bins = NULL;
}
// 1st time allocation of xhold and bins
if (dist_check) {
if (maxhold == 0) {
maxhold = atom->nmax;
xhold = memory->create_2d_double_array(maxhold,3,"neigh:xhold");
}
}
if (style != NSQ) {
if (maxbin == 0) {
maxbin = atom->nmax;
bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins");
}
}
// exclusion lists for type, group, molecule settings from neigh_modify
n = atom->ntypes;
if (nex_type == 0 && nex_group == 0 && nex_mol == 0) exclude = 0;
else exclude = 1;
if (nex_type) {
memory->destroy_2d_int_array(ex_type);
ex_type = (int **) memory->create_2d_int_array(n+1,n+1,"neigh:ex_type");
for (i = 1; i <= n; i++)
for (j = 1; j <= n; j++)
ex_type[i][j] = 0;
for (i = 0; i < nex_type; i++) {
if (ex1_type[i] <= 0 || ex1_type[i] > n ||
ex2_type[i] <= 0 || ex2_type[i] > n)
error->all("Invalid atom type in neighbor exclusion list");
ex_type[ex1_type[i]][ex2_type[i]] = 1;
ex_type[ex2_type[i]][ex1_type[i]] = 1;
}
}
if (nex_group) {
delete [] ex1_bit;
delete [] ex2_bit;
ex1_bit = new int[nex_group];
ex2_bit = new int[nex_group];
for (i = 0; i < nex_group; i++) {
ex1_bit[i] = group->bitmask[ex1_group[i]];
ex2_bit[i] = group->bitmask[ex2_group[i]];
}
}
if (nex_mol) {
delete [] ex_mol_bit;
ex_mol_bit = new int[nex_mol];
for (i = 0; i < nex_mol; i++)
ex_mol_bit[i] = group->bitmask[ex_mol_group[i]];
}
// ------------------------------------------------------------------
// pairwise lists
// test if pairwise lists need to be re-created
// no need to re-create if:
// neigh style has not changed and current requests = old requests
int same = 1;
if (style != old_style) same = 0;
if (nrequest != old_nrequest) same = 0;
else
for (i = 0; i < nrequest; i++)
if (requests[i]->identical(old_requests[i]) == 0) same = 0;
printf("SAME %d\n",same);
// if old and new are not the same, create new pairwise lists
if (!same) {
// delete old lists and create new ones
for (i = 0; i < nlist; i++) delete lists[i];
delete [] lists;
delete [] pair_build;
delete [] stencil_create;
nlist = nrequest;
lists = new NeighList*[nlist];
pair_build = new PairPtr[nlist];
stencil_create = new StencilPtr[nlist];
// create individual lists, one per request
// copy dnum setting from request to list
// pass list ptr back to requestor (except for Command class)
for (i = 0; i < nlist; i++) {
lists[i] = new NeighList(lmp,pgsize);
lists[i]->index = i;
lists[i]->dnum = requests[i]->dnum;
if (requests[i]->pair) {
Pair *pair = (Pair *) requests[i]->requestor;
pair->init_list(requests[i]->id,lists[i]);
} else if (requests[i]->fix) {
Fix *fix = (Fix *) requests[i]->requestor;
fix->init_list(requests[i]->id,lists[i]);
} else if (requests[i]->compute) {
Compute *compute = (Compute *) requests[i]->requestor;
compute->init_list(requests[i]->id,lists[i]);
}
}
// detect lists that are connected to other lists
// skip: point this list at request->otherlist, copy skip ptrs from request
// copy: point this list at request->otherlist
// half_from_full: point this list at preceeding full list
// granhistory: set preceeding list's listgranhistory to this list
// also set precedding list's ptr to FixShearHistory
// respaouter: point this list at preceeding 1/2 inner/middle lists
// pair and half and non-skip: if there is a pair full non-skip list,
// change this list to half_from_full and point at the full list
// fix/compute requests:
// kind of request = half or full
// occasional or not doesn't matter
// if non-skip pair list of same kind exists, become copy of that list
// else if request is half and non-skip pair full exists,
// become half_from_full of that list
// else do nothing and the fix/compute list will be built directly
for (i = 0; i < nlist; i++) {
if (requests[i]->skip) {
lists[i]->listskip = lists[requests[i]->otherlist];
lists[i]->iskip = requests[i]->iskip;
lists[i]->ijskip = requests[i]->ijskip;
}
if (requests[i]->copy)
lists[i]->listcopy = lists[requests[i]->otherlist];
if (requests[i]->half_from_full) lists[i]->listfull = lists[i-1];
if (requests[i]->granhistory) {
lists[i-1]->listgranhistory = lists[i];
for (int ifix = 0; ifix < modify->nfix; ifix++)
if (strcmp(modify->fix[ifix]->style,"SHEAR_HISTORY") == 0)
lists[i-1]->fix_history = (FixShearHistory *) modify->fix[ifix];
}
if (requests[i]->respaouter) {
if (requests[i-1]->respainner) {
lists[i]->respamiddle = 0;
lists[i]->listinner = lists[i-1];
} else {
lists[i]->respamiddle = 1;
lists[i]->listmiddle = lists[i-1];
lists[i]->listinner = lists[i-2];
}
}
if (requests[i]->pair && requests[i]->half && requests[i]->skip == 0) {
for (j = 0; j < nlist; j++)
- if (requests[i]->pair && requests[j]->full &&
+ if (requests[j]->pair && requests[j]->full &&
requests[j]->skip == 0) break;
if (j < nlist) {
requests[i]->half = 0;
requests[i]->half_from_full = 1;
lists[i]->listfull = lists[j];
}
}
if (requests[i]->fix || requests[i]->compute) {
for (j = 0; j < nlist; j++) {
if (requests[j]->pair && requests[j]->skip == 0 &&
requests[j]->half && requests[i]->half) break;
if (requests[j]->pair && requests[j]->skip == 0 &&
requests[j]->full && requests[i]->full) break;
}
if (j < nlist) {
requests[i]->copy = 1;
lists[i]->listcopy = lists[j];
} else {
for (j = 0; j < nlist; j++) {
if (requests[j]->pair && requests[j]->skip == 0 &&
requests[j]->full && requests[i]->half) break;
}
if (j < nlist) {
requests[i]->half = 0;
requests[i]->half_from_full = 1;
lists[i]->listfull = lists[j];
}
}
}
}
// set ptrs to pair_build and stencil_create functions for each list
// ptrs set to NULL if not set explicitly
for (i = 0; i < nlist; i++) {
choose_build(i,requests[i]);
if (style != NSQ) choose_stencil(i,requests[i]);
else stencil_create[i] = NULL;
}
// set each list's build/grow/stencil flags based on neigh request
// buildflag = 1 if its pair_build() invoked every reneighbor
// growflag = 1 if it stores atom-based arrays and pages
// stencilflag = 1 if it stores stencil arrays
for (i = 0; i < nlist; i++) {
lists[i]->buildflag = 1;
if (pair_build[i] == NULL) lists[i]->buildflag = 0;
if (requests[i]->occasional) lists[i]->buildflag = 0;
lists[i]->growflag = 1;
if (requests[i]->copy) lists[i]->growflag = 0;
lists[i]->stencilflag = 1;
if (style == NSQ) lists[i]->stencilflag = 0;
if (stencil_create[i] == NULL) lists[i]->stencilflag = 0;
}
// DEBUG: print list attributes
for (i = 0; i < nlist; i++) lists[i]->print_attributes();
// allocate atom arrays and 1st pages of lists that store them
for (i = 0; i < nlist; i++)
if (lists[i]->growflag) {
lists[i]->grow(maxlocal);
lists[i]->add_pages();
}
// setup 3 vectors of pairwise neighbor lists
// blist = lists whose pair_build() is invoked every reneighbor
// glist = lists who store atom arrays which are used every reneighbor
// slist = lists who store stencil arrays which are used every reneighbor
// blist and glist vectors are used by neighbor::build()
// slist vector is used by neighbor::setup_bins()
nblist = nglist = nslist = 0;
delete [] blist;
delete [] glist;
delete [] slist;
blist = new int[nlist];
glist = new int[nlist];
slist = new int[nlist];
for (i = 0; i < nlist; i++) {
if (lists[i]->buildflag) blist[nblist++] = i;
if (lists[i]->growflag && requests[i]->occasional == 0)
glist[nglist++] = i;
if (lists[i]->stencilflag && requests[i]->occasional == 0)
slist[nslist++] = i;
}
// DEBUG: print lists of lists
if (comm->me == 0) {
printf("Build lists = %d: ",nblist);
for (i = 0; i < nblist; i++) printf("%d ",blist[i]);
printf("\n");
printf("Grow lists = %d: ",nglist);
for (i = 0; i < nglist; i++) printf("%d ",glist[i]);
printf("\n");
printf("Stencil lists = %d: ",nslist);
for (i = 0; i < nslist; i++) printf("%d ",slist[i]);
printf("\n");
}
// reorder build vector if necessary
// relevant for lists that copy/skip/half-full from parent
// the derived lst must appear in blist after the parent list
// no occasional lists are in build vector
// swap two lists within blist when dependency is mis-ordered
// done when entire pass thru blist results in no swaps
int done = 0;
while (!done) {
done = 1;
for (i = 0; i < nblist; i++) {
NeighList *ptr = NULL;
if (lists[blist[i]]->listfull) ptr = lists[blist[i]]->listfull;
if (lists[blist[i]]->listcopy) ptr = lists[blist[i]]->listcopy;
if (lists[blist[i]]->listskip) ptr = lists[blist[i]]->listskip;
if (ptr == NULL) continue;
for (m = 0; m < nlist; m++)
if (ptr == lists[m]) break;
for (j = 0; j < nblist; j++)
if (m == blist[j]) break;
if (j < i) continue;
int tmp = blist[i];
blist[i] = blist[j];
blist[j] = tmp;
done = 0;
break;
}
}
// DEBUG: print lists of lists
if (comm->me == 0) {
printf("Build lists = %d: ",nblist);
for (i = 0; i < nblist; i++) printf("%d ",blist[i]);
printf("\n");
printf("Grow lists = %d: ",nglist);
for (i = 0; i < nglist; i++) printf("%d ",glist[i]);
printf("\n");
printf("Stencil lists = %d: ",nslist);
for (i = 0; i < nslist; i++) printf("%d ",slist[i]);
printf("\n");
}
}
// delete old requests
// copy current requests and style to old for next run
for (i = 0; i < old_nrequest; i++) delete old_requests[i];
memory->sfree(old_requests);
old_nrequest = nrequest;
old_requests = requests;
nrequest = maxrequest = 0;
requests = NULL;
old_style = style;
// ------------------------------------------------------------------
// topology lists
// 1st time allocation of topology lists
if (atom->molecular && atom->nbonds && maxbond == 0) {
if (nprocs == 1) maxbond = atom->nbonds;
else maxbond = static_cast<int> (LB_FACTOR * atom->nbonds / nprocs);
bondlist = memory->create_2d_int_array(maxbond,3,"neigh:bondlist");
}
if (atom->molecular && atom->nangles && maxangle == 0) {
if (nprocs == 1) maxangle = atom->nangles;
else maxangle = static_cast<int> (LB_FACTOR * atom->nangles / nprocs);
anglelist = memory->create_2d_int_array(maxangle,4,"neigh:anglelist");
}
if (atom->molecular && atom->ndihedrals && maxdihedral == 0) {
if (nprocs == 1) maxdihedral = atom->ndihedrals;
else maxdihedral = static_cast<int>
(LB_FACTOR * atom->ndihedrals / nprocs);
dihedrallist =
memory->create_2d_int_array(maxdihedral,5,"neigh:dihedrallist");
}
if (atom->molecular && atom->nimpropers && maximproper == 0) {
if (nprocs == 1) maximproper = atom->nimpropers;
else maximproper = static_cast<int>
(LB_FACTOR * atom->nimpropers / nprocs);
improperlist =
memory->create_2d_int_array(maximproper,5,"neigh:improperlist");
}
// set flags that determine which topology neighboring routines to use
// SHAKE sets bonds and angles negative
// bond_quartic sets bonds to 0
// delete_bonds sets all interactions negative
int bond_off = 0;
int angle_off = 0;
for (i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"shake") == 0)
bond_off = angle_off = 1;
if (force->bond && force->bond_match("quartic")) bond_off = 1;
if (atom->avec->bonds_allow) {
for (i = 0; i < atom->nlocal; i++) {
if (bond_off) break;
for (m = 0; m < atom->num_bond[i]; m++)
if (atom->bond_type[i][m] <= 0) bond_off = 1;
}
}
if (atom->avec->angles_allow) {
for (i = 0; i < atom->nlocal; i++) {
if (angle_off) break;
for (m = 0; m < atom->num_angle[i]; m++)
if (atom->angle_type[i][m] <= 0) angle_off = 1;
}
}
int dihedral_off = 0;
if (atom->avec->dihedrals_allow) {
for (i = 0; i < atom->nlocal; i++) {
if (dihedral_off) break;
for (m = 0; m < atom->num_dihedral[i]; m++)
if (atom->dihedral_type[i][m] <= 0) dihedral_off = 1;
}
}
int improper_off = 0;
if (atom->avec->impropers_allow) {
for (i = 0; i < atom->nlocal; i++) {
if (improper_off) break;
for (m = 0; m < atom->num_improper[i]; m++)
if (atom->improper_type[i][m] <= 0) improper_off = 1;
}
}
// set ptrs to topology build functions
if (bond_off) bond_build = &Neighbor::bond_partial;
else bond_build = &Neighbor::bond_all;
if (angle_off) angle_build = &Neighbor::angle_partial;
else angle_build = &Neighbor::angle_all;
if (dihedral_off) dihedral_build = &Neighbor::dihedral_partial;
else dihedral_build = &Neighbor::dihedral_all;
if (improper_off) improper_build = &Neighbor::improper_partial;
else improper_build = &Neighbor::improper_all;
// set topology neighbor list counts to 0
// in case all are turned off but potential is still defined
nbondlist = nanglelist = ndihedrallist = nimproperlist = 0;
}
/* ---------------------------------------------------------------------- */
int Neighbor::request(void *requestor)
{
if (nrequest == maxrequest) {
maxrequest += RQDELTA;
requests = (NeighRequest **)
memory->srealloc(requests,maxrequest*sizeof(NeighRequest *),
"neighbor:requests");
}
requests[nrequest] = new NeighRequest(lmp);
requests[nrequest]->requestor = requestor;
nrequest++;
return nrequest-1;
}
/* ----------------------------------------------------------------------
determine which pair_build function each neigh list needs
based on settings of neigh request
skip or copy -> single function
half_from_full, half, full, gran, respaouter -> choose by newton and tri
style NSQ options = newton off, newton on
style BIN options = newton off, newton on and not tri, newton on and tri
stlye MULTI options = same options as BIN
if none of these, ptr = NULL since pair_build is not invoked for this list
use "else if" b/c skip,copy can be set in addition to half,full,etc
------------------------------------------------------------------------- */
void Neighbor::choose_build(int index, NeighRequest *rq)
{
PairPtr pb = NULL;
if (rq->skip) pb = &Neighbor::skip_from;
else if (rq->copy) pb = &Neighbor::copy_from;
else if (rq->half_from_full) {
if (newton_pair == 0) pb = &Neighbor::half_full_no_newton;
else if (newton_pair == 1) pb = &Neighbor::half_full_newton;
} else if (rq->half) {
if (style == NSQ) {
if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton;
else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton;
} else if (style == BIN) {
if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton;
else if (triclinic == 0) pb = &Neighbor::half_bin_newton;
else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri;
} else if (style == MULTI) {
if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton;
else if (triclinic == 0) pb = &Neighbor::half_multi_newton;
else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri;
}
} else if (rq->full) {
if (style == NSQ) pb = &Neighbor::full_nsq;
else if (style == BIN) pb = &Neighbor::full_bin;
else if (style == MULTI) pb = &Neighbor::full_multi;
} else if (rq->gran) {
if (style == NSQ) {
if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton;
else if (newton_pair == 1) pb = &Neighbor::granular_nsq_newton;
} else if (style == BIN) {
if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton;
else if (triclinic == 0) pb = &Neighbor::granular_bin_newton;
else if (triclinic == 1) pb = &Neighbor::granular_bin_newton_tri;
} else if (style == MULTI)
error->all("Neighbor multi not yet enabled for granular");
} else if (rq->respaouter) {
if (style == NSQ) {
if (newton_pair == 0) pb = &Neighbor::respa_nsq_no_newton;
else if (newton_pair == 1) pb = &Neighbor::respa_nsq_newton;
} else if (style == BIN) {
if (newton_pair == 0) pb = &Neighbor::respa_bin_no_newton;
else if (triclinic == 0) pb = &Neighbor::respa_bin_newton;
else if (triclinic == 1) pb = &Neighbor::respa_bin_newton_tri;
} else if (style == MULTI)
error->all("Neighbor multi not yet enabled for rRESPA");
}
pair_build[index] = pb;
}
/* ----------------------------------------------------------------------
determine which stencil_create function each neigh list needs
based on settings of neigh request, only called if style != NSQ
skip or copy or half_from_full -> no stencil
half, gran, respaouter, full -> choose by newton and tri and dimension
if none of these, ptr = NULL since this list needs no stencils
use "else if" b/c skip,copy can be set in addition to half,full,etc
------------------------------------------------------------------------- */
void Neighbor::choose_stencil(int index, NeighRequest *rq)
{
StencilPtr sc = NULL;
if (rq->skip || rq->copy || rq->half_from_full) sc = NULL;
else if (rq->half || rq->gran || rq->respaouter) {
if (style == BIN) {
if (newton_pair == 0) {
if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_no_newton;
else if (dimension == 3) sc = &Neighbor::stencil_half_bin_3d_no_newton;
} else if (triclinic == 0) {
if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_newton;
else if (dimension == 3) sc = &Neighbor::stencil_half_bin_3d_newton;
} else if (triclinic == 1) {
if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_newton_tri;
else if (dimension == 3)
sc = &Neighbor::stencil_half_bin_3d_newton_tri;
}
} else if (style == MULTI) {
if (newton_pair == 0) {
if (dimension == 2) sc = &Neighbor::stencil_half_multi_2d_no_newton;
else if (dimension == 3)
sc = &Neighbor::stencil_half_multi_3d_no_newton;
} else if (triclinic == 0) {
if (dimension == 2) sc = &Neighbor::stencil_half_multi_2d_newton;
else if (dimension == 3) sc = &Neighbor::stencil_half_multi_3d_newton;
} else if (triclinic == 1) {
if (dimension == 2) sc = &Neighbor::stencil_half_multi_2d_newton_tri;
else if (dimension == 3)
sc = &Neighbor::stencil_half_multi_3d_newton_tri;
}
}
} else if (rq->full) {
if (style == BIN) {
if (dimension == 2) sc = &Neighbor::stencil_full_bin_2d;
else if (dimension == 3) sc = &Neighbor::stencil_full_bin_3d;
} else if (style == MULTI) {
if (dimension == 2) sc = &Neighbor::stencil_full_multi_2d;
else if (dimension == 3) sc = &Neighbor::stencil_full_multi_3d;
}
}
stencil_create[index] = sc;
}
/* ---------------------------------------------------------------------- */
int Neighbor::decide()
{
if (must_check) {
int n = update->ntimestep;
if (restart_check && n == output->next_restart) return 1;
for (int i = 0; i < fix_check; i++)
if (n == modify->fix[fixchecklist[i]]->next_reneighbor) return 1;
}
ago++;
if (ago >= delay && ago % every == 0) {
if (dist_check == 0) return 1;
else return check_distance();
} else return 0;
}
/* ---------------------------------------------------------------------- */
int Neighbor::check_distance()
{
double delx,dely,delz,rsq;
int nlocal = atom->nlocal;
double **x = atom->x;
int flag = 0;
for (int i = 0; i < nlocal; i++) {
delx = x[i][0] - xhold[i][0];
dely = x[i][1] - xhold[i][1];
delz = x[i][2] - xhold[i][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq > triggersq) flag = 1;
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flagall && ago == MAX(every,delay)) ndanger++;
return flagall;
}
/* ----------------------------------------------------------------------
build all perpetual neighbor lists every few timesteps
pairwise & topology lists are created as needed
------------------------------------------------------------------------- */
void Neighbor::build()
{
int i;
ago = 0;
ncalls++;
// store current atom positions if needed
if (dist_check) {
double **x = atom->x;
int nlocal = atom->nlocal;
if (nlocal > maxhold) {
maxhold = atom->nmax;
memory->destroy_2d_double_array(xhold);
xhold = memory->create_2d_double_array(maxhold,3,"neigh:xhold");
}
for (i = 0; i < nlocal; i++) {
xhold[i][0] = x[i][0];
xhold[i][1] = x[i][1];
xhold[i][2] = x[i][2];
}
}
// if necessary, extend atom arrays in pairwise lists
// only done for lists with growflag set which are used every reneighbor
if (atom->nlocal > maxlocal) {
maxlocal = atom->nmax;
for (i = 0; i < nglist; i++) lists[glist[i]]->grow(maxlocal);
}
// extend atom bin list if necessary
if (style != NSQ && atom->nmax > maxbin) {
maxbin = atom->nmax;
memory->sfree(bins);
bins = (int *) memory->smalloc(maxbin*sizeof(int),"bins");
}
// invoke building of pair and molecular neighbor lists
// only done for pairwise lists with buildflag set
for (i = 0; i < nblist; i++)
(this->*pair_build[blist[i]])(lists[blist[i]]);
/*
if (comm->me == 0) {
for (int m = 0; m < nlist; m++) {
int num = 0;
for (i = 0; i < lists[m]->inum; i++)
num += lists[m]->numneigh[lists[m]->ilist[i]];
printf("List %d length = %d\n",m,num);
}
}
*/
if (atom->molecular) {
if (atom->nbonds) (this->*bond_build)();
if (atom->nangles) (this->*angle_build)();
if (atom->ndihedrals) (this->*dihedral_build)();
if (atom->nimpropers) (this->*improper_build)();
}
}
/* ----------------------------------------------------------------------
build a single non-active pairwise neighbor list indexed by I
called by other classes when needed occasionally
------------------------------------------------------------------------- */
void Neighbor::build_one(int i)
{
// grow atom arrays and update stencils depending on growflag & stencilflag
if (lists[i]->growflag) lists[i]->grow(maxlocal);
if (lists[i]->stencilflag) {
lists[i]->stencil_allocate(smax,style);
(this->*stencil_create[i])(lists[i],sx,sy,sz);
}
(this->*pair_build[i])(lists[i]);
}
/* ----------------------------------------------------------------------
setup neighbor binning parameters
bin numbering in each dimension is global:
0 = 0.0 to binsize, 1 = binsize to 2*binsize, etc
nbin-1,nbin,etc = bbox-binsize to binsize, bbox to bbox+binsize, etc
-1,-2,etc = -binsize to 0.0, -2*size to -size, etc
code will work for any binsize
since next(xyz) and stencil extend as far as necessary
binsize = 1/2 of cutoff is roughly optimal
for orthogonal boxes:
a dim must be filled exactly by integer # of bins
in periodic, procs on both sides of PBC must see same bin boundary
in non-periodic, coord2bin() still assumes this by use of nbin xyz
for triclinic boxes:
tilted simulation box cannot contain integer # of bins
stencil & neigh list built differently to account for this
mbinlo = lowest global bin any of my ghost atoms could fall into
mbinhi = highest global bin any of my ghost atoms could fall into
mbin = number of bins I need in a dimension
------------------------------------------------------------------------- */
void Neighbor::setup_bins()
{
// bbox lo/hi = bounding box of entire domain
// bbox = size of bbox of entire domain
// bsubbox lo/hi = bounding box of my subdomain extended by ghost atoms
// for triclinic:
// bbox bounds all 8 corners of tilted box
// subdomain is in lamda coords
// include dimension-dependent extension via comm->cutghost
// domain->bbox() converts lamda extent to box coords and computes bbox
double bbox[3],bsubboxlo[3],bsubboxhi[3];
if (triclinic == 0) {
bboxlo = domain->boxlo;
bboxhi = domain->boxhi;
bsubboxlo[0] = domain->sublo[0] - cutghost;
bsubboxlo[1] = domain->sublo[1] - cutghost;
bsubboxlo[2] = domain->sublo[2] - cutghost;
bsubboxhi[0] = domain->subhi[0] + cutghost;
bsubboxhi[1] = domain->subhi[1] + cutghost;
bsubboxhi[2] = domain->subhi[2] + cutghost;
} else {
bboxlo = domain->boxlo_bound;
bboxhi = domain->boxhi_bound;
double lo[3],hi[3];
lo[0] = domain->sublo_lamda[0] - comm->cutghost[0];
lo[1] = domain->sublo_lamda[1] - comm->cutghost[1];
lo[2] = domain->sublo_lamda[2] - comm->cutghost[2];
hi[0] = domain->subhi_lamda[0] + comm->cutghost[0];
hi[1] = domain->subhi_lamda[1] + comm->cutghost[1];
hi[2] = domain->subhi_lamda[2] + comm->cutghost[2];
domain->bbox(lo,hi,bsubboxlo,bsubboxhi);
}
bbox[0] = bboxhi[0] - bboxlo[0];
bbox[1] = bboxhi[1] - bboxlo[1];
bbox[2] = bboxhi[2] - bboxlo[2];
// optimal bin size is roughly 1/2 the cutoff
// for BIN style, binsize = 1/2 of max neighbor cutoff
// for MULTI style, binsize = 1/2 of min neighbor cutoff
// special case of all cutoffs = 0.0, binsize = box size
double binsize_optimal;
if (style == BIN) binsize_optimal = 0.5*cutneighmax;
else binsize_optimal = 0.5*cutneighmin;
if (binsize_optimal == 0.0) binsize_optimal = bbox[0];
double binsizeinv = 1.0/binsize_optimal;
// test for too many global bins in any dimension due to huge global domain
if (bbox[0]*binsizeinv > INT_MAX || bbox[1]*binsizeinv > INT_MAX ||
bbox[2]*binsizeinv > INT_MAX)
error->all("Domain too large for neighbor bins");
// create actual bins
// always have one bin even if cutoff > bbox
// for 2d, nbinz = 1
nbinx = static_cast<int> (bbox[0]*binsizeinv);
nbiny = static_cast<int> (bbox[1]*binsizeinv);
if (dimension == 3) nbinz = static_cast<int> (bbox[2]*binsizeinv);
else nbinz = 1;
if (nbinx == 0) nbinx = 1;
if (nbiny == 0) nbiny = 1;
if (nbinz == 0) nbinz = 1;
// compute actual bin size for nbins to fit into box exactly
// error if actual bin size << cutoff, since will create a zillion bins
// this happens when nbin = 1 and box size << cutoff
// typically due to non-periodic, flat system in a particular dim
// in that extreme case, should use NSQ not BIN neighbor style
binsizex = bbox[0]/nbinx;
binsizey = bbox[1]/nbiny;
binsizez = bbox[2]/nbinz;
bininvx = 1.0 / binsizex;
bininvy = 1.0 / binsizey;
bininvz = 1.0 / binsizez;
if (binsize_optimal*bininvx > CUT2BIN_RATIO ||
binsize_optimal*bininvy > CUT2BIN_RATIO ||
binsize_optimal*bininvz > CUT2BIN_RATIO)
error->all("Cannot use neighbor bins - box size << cutoff");
// mbinlo/hi = lowest and highest global bins my ghost atoms could be in
// coord = lowest and highest values of coords for my ghost atoms
// static_cast(-1.5) = -1, so subract additional -1
// add in SMALL for round-off safety
int mbinxhi,mbinyhi,mbinzhi;
double coord;
coord = bsubboxlo[0] - SMALL*bbox[0];
mbinxlo = static_cast<int> ((coord-bboxlo[0])*bininvx);
if (coord < bboxlo[0]) mbinxlo = mbinxlo - 1;
coord = bsubboxhi[0] + SMALL*bbox[0];
mbinxhi = static_cast<int> ((coord-bboxlo[0])*bininvx);
coord = bsubboxlo[1] - SMALL*bbox[1];
mbinylo = static_cast<int> ((coord-bboxlo[1])*bininvy);
if (coord < bboxlo[1]) mbinylo = mbinylo - 1;
coord = bsubboxhi[1] + SMALL*bbox[1];
mbinyhi = static_cast<int> ((coord-bboxlo[1])*bininvy);
if (dimension == 3) {
coord = bsubboxlo[2] - SMALL*bbox[2];
mbinzlo = static_cast<int> ((coord-bboxlo[2])*bininvz);
if (coord < bboxlo[2]) mbinzlo = mbinzlo - 1;
coord = bsubboxhi[2] + SMALL*bbox[2];
mbinzhi = static_cast<int> ((coord-bboxlo[2])*bininvz);
}
// extend bins by 1 to insure stencil extent is included
// if 2d, only 1 bin in z
mbinxlo = mbinxlo - 1;
mbinxhi = mbinxhi + 1;
mbinx = mbinxhi - mbinxlo + 1;
mbinylo = mbinylo - 1;
mbinyhi = mbinyhi + 1;
mbiny = mbinyhi - mbinylo + 1;
if (dimension == 3) {
mbinzlo = mbinzlo - 1;
mbinzhi = mbinzhi + 1;
} else mbinzlo = mbinzhi = 0;
mbinz = mbinzhi - mbinzlo + 1;
// memory for bin ptrs
mbins = mbinx*mbiny*mbinz;
if (mbins > maxhead) {
maxhead = mbins;
memory->sfree(binhead);
binhead = (int *) memory->smalloc(maxhead*sizeof(int),"neigh:binhead");
}
// create stencil of bins to search over in neighbor list construction
// sx,sy,sz = max range of stencil extent
// smax =
// stencil is empty if cutneighmax = 0.0
sx = static_cast<int> (cutneighmax*bininvx);
if (sx*binsizex < cutneighmax) sx++;
sy = static_cast<int> (cutneighmax*bininvy);
if (sy*binsizey < cutneighmax) sy++;
sz = static_cast<int> (cutneighmax*bininvz);
if (sz*binsizez < cutneighmax) sz++;
if (dimension == 2) sz = 0;
smax = (2*sx+1) * (2*sy+1) * (2*sz+1);
// create stencils for pairwise neighbor lists
// only done for lists with stencilflag and buildflag set
for (int i = 0; i < nslist; i++) {
lists[slist[i]]->stencil_allocate(smax,style);
(this->*stencil_create[slist[i]])(lists[slist[i]],sx,sy,sz);
}
}
/* ----------------------------------------------------------------------
compute closest distance between central bin (0,0,0) and bin (i,j,k)
------------------------------------------------------------------------- */
double Neighbor::bin_distance(int i, int j, int k)
{
double delx,dely,delz;
if (i > 0) delx = (i-1)*binsizex;
else if (i == 0) delx = 0.0;
else delx = (i+1)*binsizex;
if (j > 0) dely = (j-1)*binsizey;
else if (j == 0) dely = 0.0;
else dely = (j+1)*binsizey;
if (k > 0) delz = (k-1)*binsizez;
else if (k == 0) delz = 0.0;
else delz = (k+1)*binsizez;
return (delx*delx + dely*dely + delz*delz);
}
/* ----------------------------------------------------------------------
set neighbor style and skin distance
------------------------------------------------------------------------- */
void Neighbor::set(int narg, char **arg)
{
if (narg != 2) error->all("Illegal neighbor command");
skin = atof(arg[0]);
if (skin < 0.0) error->all("Illegal neighbor command");
if (strcmp(arg[1],"nsq") == 0) style = NSQ;
else if (strcmp(arg[1],"bin") == 0) style = BIN;
else if (strcmp(arg[1],"multi") == 0) style = MULTI;
else error->all("Illegal neighbor command");
}
/* ----------------------------------------------------------------------
modify parameters of the pair-wise neighbor build
------------------------------------------------------------------------- */
void Neighbor::modify_params(int narg, char **arg)
{
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"every") == 0) {
if (iarg+2 > narg) error->all("Illegal neigh_modify command");
every = atoi(arg[iarg+1]);
if (every <= 0) error->all("Illegal neigh_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"delay") == 0) {
if (iarg+2 > narg) error->all("Illegal neigh_modify command");
delay = atoi(arg[iarg+1]);
if (delay < 0) error->all("Illegal neigh_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"check") == 0) {
if (iarg+2 > narg) error->all("Illegal neigh_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) dist_check = 1;
else if (strcmp(arg[iarg+1],"no") == 0) dist_check = 0;
else error->all("Illegal neigh_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"page") == 0) {
if (iarg+2 > narg) error->all("Illegal neigh_modify command");
pgsize = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"one") == 0) {
if (iarg+2 > narg) error->all("Illegal neigh_modify command");
oneatom = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"exclude") == 0) {
if (iarg+2 > narg) error->all("Illegal neigh_modify command");
if (strcmp(arg[iarg+1],"type") == 0) {
if (iarg+4 > narg) error->all("Illegal neigh_modify command");
if (nex_type == maxex_type) {
maxex_type += EXDELTA;
ex1_type = (int *) memory->srealloc(ex1_type,maxex_type*sizeof(int),
"neigh:ex1_type");
ex2_type = (int *) memory->srealloc(ex2_type,maxex_type*sizeof(int),
"neigh:ex2_type");
}
ex1_type[nex_type] = atoi(arg[iarg+2]);
ex2_type[nex_type] = atoi(arg[iarg+3]);
nex_type++;
iarg += 4;
} else if (strcmp(arg[iarg+1],"group") == 0) {
if (iarg+4 > narg) error->all("Illegal neigh_modify command");
if (nex_group == maxex_group) {
maxex_group += EXDELTA;
ex1_group =
(int *) memory->srealloc(ex1_group,maxex_group*sizeof(int),
"neigh:ex1_group");
ex2_group =
(int *) memory->srealloc(ex2_group,maxex_group*sizeof(int),
"neigh:ex2_group");
}
ex1_group[nex_group] = group->find(arg[iarg+2]);
ex2_group[nex_group] = group->find(arg[iarg+3]);
if (ex1_group[nex_group] == -1 || ex2_group[nex_group] == -1)
error->all("Invalid group ID in neigh_modify command");
nex_group++;
iarg += 4;
} else if (strcmp(arg[iarg+1],"molecule") == 0) {
if (iarg+3 > narg) error->all("Illegal neigh_modify command");
if (atom->molecular == 0) {
char *str = (char *)
"Must use molecular atom style with neigh_modify exclude molecule";
error->all(str);
}
if (nex_mol == maxex_mol) {
maxex_mol += EXDELTA;
ex_mol_group =
(int *) memory->srealloc(ex_mol_group,maxex_mol*sizeof(int),
"neigh:ex_mol_group");
}
ex_mol_group[nex_mol] = group->find(arg[iarg+2]);
if (ex_mol_group[nex_mol] == -1)
error->all("Invalid group ID in neigh_modify command");
nex_mol++;
iarg += 3;
} else if (strcmp(arg[iarg+1],"none") == 0) {
nex_type = nex_group = nex_mol = 0;
iarg += 2;
} else error->all("Illegal neigh_modify command");
} else error->all("Illegal neigh_modify command");
}
}
/* ----------------------------------------------------------------------
determine if atom j is in special list of atom i
if it is not, return 0
if it is and special flag is 0 (both coeffs are 0.0), return -1
if it is and special flag is 1 (both coeffs are 1.0), return 0
if it is and special flag is 2 (otherwise), return 1,2,3
for which neighbor it is (and which coeff it maps to)
------------------------------------------------------------------------- */
int Neighbor::find_special(int i, int j)
{
int *list = atom->special[i];
int n1 = atom->nspecial[i][0];
int n2 = atom->nspecial[i][1];
int n3 = atom->nspecial[i][2];
int tag = atom->tag[j];
for (int i = 0; i < n3; i++) {
if (list[i] == tag) {
if (i < n1) {
if (special_flag[1] == 0) return -1;
else if (special_flag[1] == 1) return 0;
else return 1;
} else if (i < n2) {
if (special_flag[2] == 0) return -1;
else if (special_flag[2] == 1) return 0;
else return 2;
} else {
if (special_flag[3] == 0) return -1;
else if (special_flag[3] == 1) return 0;
else return 3;
}
}
}
return 0;
}
/* ----------------------------------------------------------------------
bin owned and ghost atoms
------------------------------------------------------------------------- */
void Neighbor::bin_atoms()
{
int i,ibin,nlocal,nall;
double **x;
nlocal = atom->nlocal;
nall = atom->nlocal + atom->nghost;
x = atom->x;
for (i = 0; i < mbins; i++) binhead[i] = -1;
// bin in reverse order so linked list will be in forward order
// also puts ghost atoms at end of list, which is necessary
for (i = nall-1; i >= 0; i--) {
ibin = coord2bin(x[i]);
bins[i] = binhead[ibin];
binhead[ibin] = i;
}
}
/* ----------------------------------------------------------------------
convert atom coords into local bin #
for orthogonal, only ghost atoms will have coord >= bboxhi or coord < bboxlo
take special care to insure ghosts are put in correct bins
this is necessary so that both procs on either side of PBC
treat a pair of atoms straddling the PBC in a consistent way
for triclinic, doesn't matter since stencil & neigh list built differently
------------------------------------------------------------------------- */
int Neighbor::coord2bin(double *x)
{
int ix,iy,iz;
if (x[0] >= bboxhi[0])
ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx) + nbinx - mbinxlo;
else if (x[0] >= bboxlo[0])
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - mbinxlo;
else
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - mbinxlo - 1;
if (x[1] >= bboxhi[1])
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy) + nbiny - mbinylo;
else if (x[1] >= bboxlo[1])
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - mbinylo;
else
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - mbinylo - 1;
if (x[2] >= bboxhi[2])
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz) + nbinz - mbinzlo;
else if (x[2] >= bboxlo[2])
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - mbinzlo;
else
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - mbinzlo - 1;
return (iz*mbiny*mbinx + iy*mbinx + ix + 1);
}
/* ----------------------------------------------------------------------
test if atom pair i,j is excluded from neighbor list
due to type, group, molecule settings from neigh_modify command
return 1 if should be excluded, 0 if included
------------------------------------------------------------------------- */
int Neighbor::exclusion(int i, int j, int itype, int jtype,
int *mask, int *molecule)
{
int m;
if (nex_type && ex_type[itype][jtype]) return 1;
if (nex_group) {
for (m = 0; m < nex_group; m++) {
if (mask[i] & ex1_bit[m] && mask[j] & ex2_bit[m]) return 1;
if (mask[i] & ex2_bit[m] && mask[j] & ex1_bit[m]) return 1;
}
}
if (nex_mol) {
for (m = 0; m < nex_mol; m++)
if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] &&
molecule[i] == molecule[j]) return 1;
}
return 0;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
------------------------------------------------------------------------- */
double Neighbor::memory_usage()
{
double bytes = 0.0;
bytes += maxhold*3 * sizeof(double);
if (style != NSQ) {
bytes += maxbin * sizeof(int);
bytes += maxhead * sizeof(int);
}
for (int i = 0; i < nlist; i++) bytes += lists[i]->memory_usage();
bytes += maxbond*3 * sizeof(int);
bytes += maxangle*4 * sizeof(int);
bytes += maxdihedral*5 * sizeof(int);
bytes += maximproper*5 * sizeof(int);
return bytes;
}