From 10cb6b72e9c02085ee03399ca10b75e7fb52db9d Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 3 Oct 2007 16:26:20 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@938 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_attribute_atom.cpp | 191 ++++++++++++++++++ src/compute_attribute_atom.h | 36 ++++ src/compute_sum_atom.cpp | 158 +++++++++++++++ src/compute_sum_atom.h | 37 ++++ src/fix_ave_atom.cpp | 267 +++++++++++++++++++++++++ src/fix_ave_atom.h | 48 +++++ src/neigh_derive.cpp | 244 +++++++++++++++++++++++ src/neigh_half_multi.cpp | 347 +++++++++++++++++++++++++++++++++ src/neigh_half_nsq.cpp | 192 ++++++++++++++++++ src/neigh_list.cpp | 257 ++++++++++++++++++++++++ src/neigh_list.h | 88 +++++++++ src/neigh_request.cpp | 164 ++++++++++++++++ src/neigh_request.h | 83 ++++++++ src/pair_hybrid_overlay.cpp | 140 +++++++++++++ src/pair_hybrid_overlay.h | 32 +++ 15 files changed, 2284 insertions(+) create mode 100644 src/compute_attribute_atom.cpp create mode 100644 src/compute_attribute_atom.h create mode 100644 src/compute_sum_atom.cpp create mode 100644 src/compute_sum_atom.h create mode 100644 src/fix_ave_atom.cpp create mode 100644 src/fix_ave_atom.h create mode 100644 src/neigh_derive.cpp create mode 100644 src/neigh_half_multi.cpp create mode 100644 src/neigh_half_nsq.cpp create mode 100644 src/neigh_list.cpp create mode 100644 src/neigh_list.h create mode 100644 src/neigh_request.cpp create mode 100644 src/neigh_request.h create mode 100644 src/pair_hybrid_overlay.cpp create mode 100644 src/pair_hybrid_overlay.h diff --git a/src/compute_attribute_atom.cpp b/src/compute_attribute_atom.cpp new file mode 100644 index 0000000000..16eeaad34d --- /dev/null +++ b/src/compute_attribute_atom.cpp @@ -0,0 +1,191 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_attribute_atom.h" +#include "atom.h" +#include "domain.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{X,Y,Z,XU,YU,ZU,VX,VY,VZ,FX,FY,FZ,XYZ,V,F}; + +/* --------------------------------------------------------------------- */ + +ComputeAttributeAtom::ComputeAttributeAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all("Illegal compute ke/atom command"); + + peratom_flag = 1; + size_peratom = 0; + allocate = 1; + + if (strcmp(arg[3],"x") == 0) which = X; + else if (strcmp(arg[3],"y") == 0) which = Y; + else if (strcmp(arg[3],"z") == 0) which = Z; + else if (strcmp(arg[3],"xu") == 0) which = XU; + else if (strcmp(arg[3],"yu") == 0) which = YU; + else if (strcmp(arg[3],"zu") == 0) which = ZU; + else if (strcmp(arg[3],"vx") == 0) which = VX; + else if (strcmp(arg[3],"vy") == 0) which = VY; + else if (strcmp(arg[3],"vz") == 0) which = VZ; + else if (strcmp(arg[3],"fx") == 0) which = FX; + else if (strcmp(arg[3],"fy") == 0) which = FY; + else if (strcmp(arg[3],"fz") == 0) which = FZ; + + else if (strcmp(arg[3],"xyz") == 0) { + which = XYZ; + size_peratom = 3; + allocate = 0; + } else if (strcmp(arg[3],"v") == 0) { + which = V; + size_peratom = 3; + allocate = 0; + } else if (strcmp(arg[3],"f") == 0) { + which = F; + size_peratom = 3; + allocate = 0; + } else error->all("Illegal compute attribute/atom command"); + + nmax = 0; + s_attribute = NULL; + v_attribute = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeAttributeAtom::~ComputeAttributeAtom() +{ + if (allocate) { + memory->sfree(s_attribute); + memory->destroy_2d_double_array(v_attribute); + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeAttributeAtom::compute_peratom() +{ + // grow attribute array if necessary + + if (allocate && atom->nlocal > nmax) { + if (size_peratom == 0) { + memory->sfree(s_attribute); + nmax = atom->nmax; + s_attribute = (double *) + memory->smalloc(nmax*sizeof(double), + "compute/attribute/atom:s_attribute"); + } else { + memory->destroy_2d_double_array(v_attribute); + nmax = atom->nmax; + v_attribute = + memory->create_2d_double_array(nmax,size_peratom, + "compute/attribute/atom:v_attribute"); + } + } + + // fill attribute vector with appropriate atom value + // or simply set pointer to exisitng atom vector + + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + int *mask = atom->mask; + int *image = atom->image; + int nlocal = atom->nlocal; + + if (which == X) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_attribute[i] = x[i][0]; + else s_attribute[i] = 0.0; + } else if (which == Y) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_attribute[i] = x[i][1]; + else s_attribute[i] = 0.0; + } else if (which == Z) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_attribute[i] = x[i][2]; + else s_attribute[i] = 0.0; + + } else if (which == XU) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + s_attribute[i] = x[i][0] + ((image[i] & 1023) - 512) * xprd; + else s_attribute[i] = 0.0; + } else if (which == YU) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + s_attribute[i] = x[i][1] + ((image[i] >> 10 & 1023) - 512) * yprd; + else s_attribute[i] = 0.0; + } else if (which == ZU) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + s_attribute[i] = x[i][2] + ((image[i] >> 20) - 512) * zprd; + else s_attribute[i] = 0.0; + + } else if (which == VX) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_attribute[i] = v[i][0]; + else s_attribute[i] = 0.0; + } else if (which == VY) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_attribute[i] = v[i][1]; + else s_attribute[i] = 0.0; + } else if (which == VZ) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_attribute[i] = v[i][2]; + else s_attribute[i] = 0.0; + + } else if (which == FX) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_attribute[i] = f[i][0]; + else s_attribute[i] = 0.0; + } else if (which == FY) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_attribute[i] = f[i][1]; + else s_attribute[i] = 0.0; + } else if (which == FZ) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_attribute[i] = f[i][2]; + else s_attribute[i] = 0.0; + + } else if (which == XYZ) v_attribute = x; + else if (which == V) v_attribute = v; + else if (which == F) v_attribute = f; + + // set appropriate compute ptr to local array + + if (size_peratom == 0) scalar_atom = s_attribute; + else vector_atom = v_attribute; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +int ComputeAttributeAtom::memory_usage() +{ + int bytes = 0; + if (allocate) { + if (size_peratom == 0) bytes = nmax * sizeof(double); + else bytes = size_peratom * nmax * sizeof(double); + } + return bytes; +} diff --git a/src/compute_attribute_atom.h b/src/compute_attribute_atom.h new file mode 100644 index 0000000000..f83819e65e --- /dev/null +++ b/src/compute_attribute_atom.h @@ -0,0 +1,36 @@ +/* ---------------------------------------------------------------------- + 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_ATTRIBUTE_ATOM_H +#define COMPUTE_ATTRIBUTE_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeAttributeAtom : public Compute { + public: + ComputeAttributeAtom(class LAMMPS *, int, char **); + ~ComputeAttributeAtom(); + void init() {} + void compute_peratom(); + int memory_usage(); + + private: + int which,allocate,nmax; + double *s_attribute,**v_attribute; +}; + +} + +#endif diff --git a/src/compute_sum_atom.cpp b/src/compute_sum_atom.cpp new file mode 100644 index 0000000000..d9b3b29711 --- /dev/null +++ b/src/compute_sum_atom.cpp @@ -0,0 +1,158 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_sum_atom.h" +#include "atom.h" +#include "modify.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeSumAtom::ComputeSumAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 5) error->all("Illegal compute sum/atom command"); + + // store pre-compute IDs + + npre = narg - 3; + id_pre = new char*[npre]; + for (int i = 0; i < npre; i++) { + int iarg = i + 3; + int n = strlen(arg[iarg]) + 1; + id_pre[i] = new char[n]; + strcpy(id_pre[i],arg[iarg]); + } + + compute = new Compute*[npre]; + + // check consistency of set of pre-computes for scalar & vector output + + peratom_flag = 1; + int icompute = modify->find_compute(id_pre[0]); + size_peratom = modify->compute[icompute]->size_peratom; + + for (int i = 1; i < npre; i++) { + icompute = modify->find_compute(id_pre[i]); + if (icompute < 0) + error->all("Could not find compute sum/atom pre-compute ID"); + if (modify->compute[icompute]->peratom_flag == 0) + error->all("Compute sum/atom compute does not compute vector per atom"); + if (modify->compute[icompute]->size_peratom != size_peratom) + error->all("Inconsistent sizes of compute sum/atom compute vectors"); + } + + nmax = 0; + s_value = NULL; + v_value = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeSumAtom::~ComputeSumAtom() +{ + delete [] compute; + memory->sfree(s_value); + memory->destroy_2d_double_array(v_value); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSumAtom::init() +{ + // set ptrs to Computes used as pre-computes by this compute + + for (int i = 0; i < npre; i++) { + int icompute = modify->find_compute(id_pre[i]); + if (icompute < 0) + error->all("Could not find compute sum/atom pre-compute ID"); + compute[i] = modify->compute[icompute]; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeSumAtom::compute_peratom() +{ + int i,j,m; + + // grow sum array if necessary + + if (atom->nlocal > nmax) { + nmax = atom->nmax; + if (size_peratom == 0) { + memory->sfree(s_value); + s_value = (double *) + memory->smalloc(nmax*sizeof(double),"compute/sum/atom:s_value"); + scalar_atom = s_value; + } else { + memory->destroy_2d_double_array(v_value); + v_value = memory->create_2d_double_array(nmax,size_peratom, + "compute/sum/atom:v_value"); + vector_atom = v_value; + } + } + + // sum over pre-computes + // pre-computes of the pre-computes are not invoked + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (size_peratom == 0) { + double *scalar = compute[0]->scalar_atom; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_value[i] = scalar[i]; + else s_value[i] = 0.0; + + for (m = 1; m < npre; m++) { + scalar = compute[m]->scalar_atom; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) s_value[i] += scalar[i]; + } + + } else { + double **vector = compute[0]->vector_atom; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + for (j = 0; j < size_peratom; j++) + v_value[i][j] = vector[i][j]; + else + for (j = 0; j < size_peratom; j++) + v_value[i][j] = 0.0; + + for (m = 1; m < npre; m++) { + vector = compute[m]->vector_atom; + for (j = 0; j < size_peratom; j++) + if (mask[i] & groupbit) v_value[i][j] += vector[i][j]; + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +int ComputeSumAtom::memory_usage() +{ + int bytes = 0; + if (size_peratom == 0) bytes = nmax * sizeof(double); + else bytes = nmax*size_peratom * sizeof(double); + return bytes; +} diff --git a/src/compute_sum_atom.h b/src/compute_sum_atom.h new file mode 100644 index 0000000000..4c4d4dd65c --- /dev/null +++ b/src/compute_sum_atom.h @@ -0,0 +1,37 @@ +/* ---------------------------------------------------------------------- + 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_SUM_ATOM_H +#define COMPUTE_SUM_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeSumAtom : public Compute { + public: + ComputeSumAtom(class LAMMPS *, int, char **); + ~ComputeSumAtom(); + void init(); + void compute_peratom(); + int memory_usage(); + + private: + int nmax; + class Compute **compute; + double *s_value,**v_value; +}; + +} + +#endif diff --git a/src/fix_ave_atom.cpp b/src/fix_ave_atom.cpp new file mode 100644 index 0000000000..a15e3c6ffd --- /dev/null +++ b/src/fix_ave_atom.cpp @@ -0,0 +1,267 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "string.h" +#include "fix_ave_atom.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "compute.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 7) error->all("Illegal fix ave/atom command"); + + nevery = atoi(arg[3]); + nrepeat = atoi(arg[4]); + peratom_freq = atoi(arg[5]); + + int n = strlen(arg[6]) + 1; + id_compute = new char[n]; + strcpy(id_compute,arg[6]); + + // setup and error check + + if (nevery <= 0) error->all("Illegal fix ave/atom command"); + if (peratom_freq < nevery || peratom_freq % nevery || + (nrepeat-1)*nevery >= peratom_freq) + error->all("Illegal fix ave/atom command"); + + int icompute = modify->find_compute(id_compute); + if (icompute < 0) error->all("Compute ID for fix ave/atom does not exist"); + + if (modify->compute[icompute]->peratom_flag == 0) + error->all("Fix ave/atom compute does not calculate per-atom info"); + + if (modify->compute[icompute]->pressflag) pressure_every = nevery; + + peratom_flag = 1; + + // setup list of computes to call, including pre-computes + + ncompute = 1 + modify->compute[icompute]->npre; + compute = new Compute*[ncompute]; + + // perform initial allocation of atom-based array + // register with Atom class + + size_peratom = modify->compute[icompute]->size_peratom; + scalar = NULL; + vector = NULL; + grow_arrays(atom->nmax); + atom->add_callback(0); + + // zero the array since dump may access it on timestep 0 + + int nlocal = atom->nlocal; + if (size_peratom == 0) + for (int i = 0; i < nlocal; i++) scalar[i] = 0.0; + else + for (int i = 0; i < nlocal; i++) + for (int m = 0; m < size_peratom; m++) + vector[i][m] = 0.0; + + // nvalid = next step on which end_of_step does something + + irepeat = 0; + nvalid = (update->ntimestep/peratom_freq)*peratom_freq + peratom_freq; + nvalid -= (nrepeat-1)*nevery; + if (nvalid <= update->ntimestep) + error->all("Fix ave/atom cannot be started on this timestep"); +} + +/* ---------------------------------------------------------------------- */ + +FixAveAtom::~FixAveAtom() +{ + // unregister callback to this fix from Atom class + + atom->delete_callback(id,0); + + delete [] id_compute; + delete [] compute; + memory->sfree(scalar); + memory->destroy_2d_double_array(vector); +} + +/* ---------------------------------------------------------------------- */ + +int FixAveAtom::setmask() +{ + int mask = 0; + mask |= END_OF_STEP; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixAveAtom::init() +{ + // set ptrs to one or more computes called each end-of-step + + int icompute = modify->find_compute(id_compute); + if (icompute < 0) error->all("Compute ID for fix ave/atom does not exist"); + + ncompute = 0; + if (modify->compute[icompute]->npre) + for (int i = 0; i < modify->compute[icompute]->npre; i++) { + int ic = modify->find_compute(modify->compute[icompute]->id_pre[i]); + if (ic < 0) + error->all("Precompute ID for fix ave/atom does not exist"); + compute[ncompute++] = modify->compute[ic]; + } + + compute[ncompute++] = modify->compute[icompute]; +} + +/* ---------------------------------------------------------------------- */ + +void FixAveAtom::end_of_step() +{ + int i,m; + + // skip if not step which requires doing something + + if (update->ntimestep != nvalid) return; + + // zero if first step + + int nlocal = atom->nlocal; + + if (irepeat == 0) { + if (size_peratom == 0) + for (i = 0; i < nlocal; i++) scalar[i] = 0.0; + else + for (i = 0; i < nlocal; i++) + for (m = 0; m < size_peratom; m++) + vector[i][m] = 0.0; + } + + // accumulate results of compute to local copy + + for (i = 0; i < ncompute; i++) compute[i]->compute_peratom(); + + int *mask = atom->mask; + + if (size_peratom == 0) { + double *compute_scalar = compute[ncompute-1]->scalar_atom; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) scalar[i] += compute_scalar[i]; + } else { + double **compute_vector = compute[ncompute-1]->vector_atom; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + for (m = 0; m < size_peratom; m++) + vector[i][m] += compute_vector[i][m]; + } + + irepeat++; + nvalid += nevery; + + // divide by nrepeat if final step + // reset irepeat and nvalid + + if (irepeat == nrepeat) { + double repeat = nrepeat; + + if (size_peratom == 0) + for (i = 0; i < nlocal; i++) + scalar[i] /= repeat; + else + for (i = 0; i < nlocal; i++) + for (m = 0; m < size_peratom; m++) + vector[i][m] /= repeat; + + irepeat = 0; + nvalid = update->ntimestep+peratom_freq - (nrepeat-1)*nevery; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +int FixAveAtom::memory_usage() +{ + int bytes; + if (size_peratom == 0) bytes = atom->nmax * sizeof(double); + else bytes = atom->nmax*size_peratom * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + allocate atom-based array +------------------------------------------------------------------------- */ + +void FixAveAtom::grow_arrays(int nmax) +{ + if (size_peratom == 0) { + scalar = (double *) memory->srealloc(scalar,nmax,"fix_ave/atom:scalar"); + scalar_atom = scalar; + } else { + vector = memory->grow_2d_double_array(vector,nmax,size_peratom, + "fix_ave/atom:vector"); + vector_atom = vector; + } +} + +/* ---------------------------------------------------------------------- + copy values within local atom-based array +------------------------------------------------------------------------- */ + +void FixAveAtom::copy_arrays(int i, int j) +{ + if (size_peratom == 0) + scalar[j] = scalar[i]; + else + for (int m = 0; m <= size_peratom; m++) + vector[j][m] = vector[i][m]; +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based array for exchange with another proc +------------------------------------------------------------------------- */ + +int FixAveAtom::pack_exchange(int i, double *buf) +{ + if (size_peratom == 0) { + buf[0] = scalar[i]; + return 1; + } + + for (int m = 0; m <= size_peratom; m++) buf[m] = vector[i][m]; + return size_peratom; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based array from exchange with another proc +------------------------------------------------------------------------- */ + +int FixAveAtom::unpack_exchange(int nlocal, double *buf) +{ + if (size_peratom == 0) { + scalar[nlocal] = buf[0]; + return 1; + } + + for (int m = 0; m <= size_peratom; m++) vector[nlocal][m] = buf[m]; + return size_peratom; +} diff --git a/src/fix_ave_atom.h b/src/fix_ave_atom.h new file mode 100644 index 0000000000..6b806e1652 --- /dev/null +++ b/src/fix_ave_atom.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + 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 FIX_AVE_ATOM_H +#define FIX_AVE_ATOM_H + +#include "stdio.h" +#include "fix.h" + +namespace LAMMPS_NS { + +class FixAveAtom : public Fix { + public: + FixAveAtom(class LAMMPS *, int, char **); + ~FixAveAtom(); + int setmask(); + void init(); + void end_of_step(); + + int memory_usage(); + void grow_arrays(int); + void copy_arrays(int, int); + int pack_exchange(int, double *); + int unpack_exchange(int, double *); + + private: + int nrepeat,irepeat,nvalid; + char *id_compute; + int ncompute; + class Compute **compute; + + double *scalar; + double **vector; +}; + +} + +#endif diff --git a/src/neigh_derive.cpp b/src/neigh_derive.cpp new file mode 100644 index 0000000000..d07176e35a --- /dev/null +++ b/src/neigh_derive.cpp @@ -0,0 +1,244 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neigh_list.h" +#include "atom.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + build half list from full list + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) + works if full list is a skip list +------------------------------------------------------------------------- */ + +void Neighbor::half_full_no_newton(NeighList *list) +{ + int i,j,ii,jj,n,jnum; + int *neighptr,*jlist; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int *ilist_full = list->listfull->ilist; + int *numneigh_full = list->listfull->numneigh; + int **firstneigh_full = list->listfull->firstneigh; + int inum_full = list->listfull->inum; + + int inum = 0; + int npage = 0; + int npnt = 0; + + // loop over atoms in full list + + for (ii = 0; ii < inum_full; ii++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + // loop over full neighbor list + + i = ilist_full[ii]; + jlist = firstneigh_full[i]; + jnum = numneigh_full[i]; + + for (j = 0; jj < jnum; jj++) { + j = jlist[jj]; + if (j > i) neighptr[n++] = j; + } + + ilist[inum] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + inum++; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = inum; +} + +/* ---------------------------------------------------------------------- + build half list from full list + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) + works if full list is a skip list +------------------------------------------------------------------------- */ + +void Neighbor::half_full_newton(NeighList *list) +{ + int i,j,ii,jj,n,jnum; + int *neighptr,*jlist; + double xtmp,ytmp,ztmp; + + double **x = atom->x; + int nlocal = atom->nlocal; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int *ilist_full = list->listfull->ilist; + int *numneigh_full = list->listfull->numneigh; + int **firstneigh_full = list->listfull->firstneigh; + int inum_full = list->listfull->inum; + + int inum = 0; + int npage = 0; + int npnt = 0; + + // loop over atoms in full list + + for (ii = 0; ii < inum_full; ii++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + i = ilist_full[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over full neighbor list + + jlist = firstneigh_full[i]; + jnum = numneigh_full[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + if (j < nlocal) { + if (i > j) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + neighptr[n++] = j; + } + + ilist[inum] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + inum++; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = inum; +} + +/* ---------------------------------------------------------------------- + build skip list for subset of types from parent list + iskip and ijskip flag which atom types and type pairs to skip +------------------------------------------------------------------------- */ + +void Neighbor::skip_from(NeighList *list) +{ + int i,j,ii,jj,n,itype,jnum,joriginal; + int *neighptr,*jlist; + + int *type = atom->type; + int nall = atom->nlocal + atom->nghost; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int *ilist_skip = list->listskip->ilist; + int *numneigh_skip = list->listskip->numneigh; + int **firstneigh_skip = list->listskip->firstneigh; + int inum_skip = list->listskip->inum; + + int *iskip = list->iskip; + int **ijskip = list->ijskip; + + int inum = 0; + int npage = 0; + int npnt = 0; + + // loop over atoms in other list + // skip I atom entirely if iskip is set for type[I] + // skip I,J pair if ijskip is set for type[I],type[J] + + for (ii = 0; ii < inum_skip; ii++) { + i = ilist_skip[ii]; + itype = type[i]; + if (iskip[type[i]]) continue; + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + // loop over full neighbor list + + jlist = firstneigh_skip[i]; + jnum = numneigh_skip[i]; + + for (jj = 0; jj < jnum; jj++) { + j = joriginal = jlist[jj]; + if (j >= nall) j %= nall; + if (ijskip[itype][type[j]]) continue; + neighptr[n++] = joriginal; + } + + ilist[inum] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + inum++; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = inum; +} + +/* ---------------------------------------------------------------------- + create list which is simply a copy of another list +------------------------------------------------------------------------- */ + +void Neighbor::copy_from(NeighList *list) +{ + NeighList *listcopy = list->listcopy; + + list->inum = listcopy->inum; + list->ilist = listcopy->ilist; + list->numneigh = listcopy->numneigh; + list->firstneigh = listcopy->firstneigh; + list->firstdouble = listcopy->firstdouble; + list->pages = listcopy->pages; + list->dpages = listcopy->dpages; +} diff --git a/src/neigh_half_multi.cpp b/src/neigh_half_multi.cpp new file mode 100644 index 0000000000..8ecad243a2 --- /dev/null +++ b/src/neigh_half_multi.cpp @@ -0,0 +1,347 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neigh_list.h" +#include "atom.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + binned neighbor list construction with partial Newton's 3rd law + each owned atom i checks own bin and other bins in stencil + multi-type stencil is itype dependent and is distance checked + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::half_multi_no_newton(NeighList *list) +{ + int i,j,k,n,itype,jtype,ibin,which,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + double *cutsq,*distsq; + + // bin local & ghost atoms + + bin_atoms(); + + // loop over each atom, storing neighbors + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int *nstencil_multi = list->nstencil_multi; + int **stencil_multi = list->stencil_multi; + double **distsq_multi = list->distsq_multi; + + int inum = 0; + int npage = 0; + int npnt = 0; + + for (i = 0; i < nlocal; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in other bins in stencil including self + // only store pair if i < j + // skip if i,j neighbor cutoff is less than bin distance + // stores own/own pairs only once + // stores own/ghost pairs on both procs + + ibin = coord2bin(x[i]); + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) 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 <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(i,j); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + } + + ilist[inum] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + inum++; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = inum; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with full Newton's 3rd law + each owned atom i checks its own bin and other bins in Newton stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::half_multi_newton(NeighList *list) +{ + int i,j,k,n,itype,jtype,ibin,which,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + double *cutsq,*distsq; + + // bin local & ghost atoms + + bin_atoms(); + + // loop over each atom, storing neighbors + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int *nstencil_multi = list->nstencil_multi; + int **stencil_multi = list->stencil_multi; + double **distsq_multi = list->distsq_multi; + + int inum = 0; + int npage = 0; + int npnt = 0; + + for (i = 0; i < nlocal; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over rest of atoms in i's bin, ghosts are at end of linked list + // if j is owned atom, store it, since j is beyond i in linked list + // if j is ghost, only store if j coords are "above and to the right" of i + + for (j = bins[i]; j >= 0; j = bins[j]) { + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) 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 <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(i,j); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + + // loop over all atoms in other bins in stencil, store every pair + // skip if i,j neighbor cutoff is less than bin distance + + ibin = coord2bin(x[i]); + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) 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 <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(i,j); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + } + + ilist[inum] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + inum++; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = inum; +} + +/* ---------------------------------------------------------------------- + binned neighbor list construction with Newton's 3rd law for triclinic + each owned atom i checks its own bin and other bins in triclinic stencil + multi-type stencil is itype dependent and is distance checked + every pair stored exactly once by some processor +------------------------------------------------------------------------- */ + +void Neighbor::half_multi_newton_tri(NeighList *list) +{ + int i,j,k,n,itype,jtype,ibin,which,ns; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr,*s; + double *cutsq,*distsq; + + // bin local & ghost atoms + + bin_atoms(); + + // loop over each atom, storing neighbors + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int *nstencil_multi = list->nstencil_multi; + int **stencil_multi = list->stencil_multi; + double **distsq_multi = list->distsq_multi; + + int inum = 0; + int npage = 0; + int npnt = 0; + + for (i = 0; i < nlocal; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in bins, including self, in stencil + // skip if i,j neighbor cutoff is less than bin distance + // bins below self are excluded from stencil + // pairs for atoms j below i are excluded + + ibin = coord2bin(x[i]); + s = stencil_multi[itype]; + distsq = distsq_multi[itype]; + cutsq = cutneighsq[itype]; + ns = nstencil_multi[itype]; + for (k = 0; k < ns; k++) { + for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) { + jtype = type[j]; + if (cutsq[jtype] < distsq[k]) continue; + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue; + + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) 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 <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(i,j); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + } + + ilist[inum] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + inum++; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = inum; +} diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp new file mode 100644 index 0000000000..f3f75698f8 --- /dev/null +++ b/src/neigh_half_nsq.cpp @@ -0,0 +1,192 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "neighbor.h" +#include "neigh_list.h" +#include "atom.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- + N^2 / 2 search for neighbor pairs with partial Newton's 3rd law + pair stored once if i,j are both owned and i < j + pair stored by me if j is ghost (also stored by proc owning j) +------------------------------------------------------------------------- */ + +void Neighbor::half_nsq_no_newton(NeighList *list) +{ + int i,j,n,itype,jtype,which; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + int inum = 0; + int npage = 0; + int npnt = 0; + + for (i = 0; i < nlocal; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over remaining atoms, owned and ghost + + for (j = i+1; j < nall; j++) { + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) 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 <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(i,j); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + + ilist[inum] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + inum++; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = inum; +} + +/* ---------------------------------------------------------------------- + N^2 / 2 search for neighbor pairs with full Newton's 3rd law + every pair stored exactly once by some processor + decision on ghost atoms based on itag,jtag tests +------------------------------------------------------------------------- */ + +void Neighbor::half_nsq_newton(NeighList *list) +{ + int i,j,n,itype,jtype,itag,jtag,which; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *neighptr; + + double **x = atom->x; + int *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = atom->nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int nstencil = list->nstencil; + int *stencil = list->stencil; + + int inum = 0; + int npage = 0; + int npnt = 0; + + for (i = 0; i < nlocal; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itag = tag[i]; + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over remaining atoms, owned and ghost + // itag = jtag is possible for long cutoffs that include images of self + + for (j = i+1; j < nall; j++) { + if (j >= nlocal) { + jtag = tag[j]; + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + else if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) + continue; + } + } + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) 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 <= cutneighsq[itype][jtype]) { + if (molecular) which = find_special(i,j); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = which*nall + j; + } + } + + ilist[inum] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + inum++; + npnt += n; + if (npnt >= pgsize) + error->one("Neighbor list overflow, boost neigh_modify one or page"); + } + + list->inum = inum; +} diff --git a/src/neigh_list.cpp b/src/neigh_list.cpp new file mode 100644 index 0000000000..10bc41429c --- /dev/null +++ b/src/neigh_list.cpp @@ -0,0 +1,257 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "neigh_list.h" +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define PGDELTA 1 + +enum{NSQ,BIN,MULTI}; // also in neighbor.cpp + +/* ---------------------------------------------------------------------- */ + +NeighList::NeighList(LAMMPS *lmp, int size) : Pointers(lmp) +{ + maxlocal = 0; + pgsize = size; + + ilist = NULL; + numneigh = NULL; + firstneigh = NULL; + firstdouble = NULL; + + maxpage = 0; + pages = NULL; + dpages = NULL; + dnum = 0; + + skip = 0; + iskip = NULL; + ijskip = NULL; + + listgranhistory = NULL; + fix_history = NULL; + + respamiddle = 0; + listinner = NULL; + listmiddle = NULL; + listfull = NULL; + listcopy = NULL; + listskip = NULL; + + maxstencil = 0; + stencil = NULL; + + maxstencil_multi = 0; + nstencil_multi = NULL; + stencil_multi = NULL; + distsq_multi = NULL; +} + +/* ---------------------------------------------------------------------- */ + +NeighList::~NeighList() +{ + if (!listcopy) { + memory->sfree(ilist); + memory->sfree(numneigh); + memory->sfree(firstneigh); + memory->sfree(firstdouble); + + for (int i = 0; i < maxpage; i++) memory->sfree(pages[i]); + memory->sfree(pages); + if (dnum) { + for (int i = 0; i < maxpage; i++) memory->sfree(dpages[i]); + memory->sfree(dpages); + } + } + + if (maxstencil) memory->sfree(stencil); + if (maxstencil_multi) { + for (int i = 1; i <= atom->ntypes; i++) { + memory->sfree(stencil_multi[i]); + memory->sfree(distsq_multi[i]); + } + delete [] nstencil_multi; + delete [] stencil_multi; + delete [] distsq_multi; + } +} + +/* ---------------------------------------------------------------------- + grow atom arrays to allow for nmax atoms + triggered by more atoms on a processor +------------------------------------------------------------------------- */ + +void NeighList::grow(int nmax) +{ + // skip if grow not needed (don't think this should happen) + + if (nmax <= maxlocal) return; + maxlocal = nmax; + + memory->sfree(ilist); + memory->sfree(numneigh); + memory->sfree(firstneigh); + ilist = (int *) + memory->smalloc(maxlocal*sizeof(int),"neighlist:ilist"); + numneigh = (int *) + memory->smalloc(maxlocal*sizeof(int),"neighlist:numneigh"); + firstneigh = (int **) + memory->smalloc(maxlocal*sizeof(int *),"neighlist:firstneigh"); + + if (dnum) + firstdouble = (double **) + memory->smalloc(maxlocal*sizeof(double *),"neighlist:firstdouble"); +} + +/* ---------------------------------------------------------------------- + insure stencils are large enough for smax bins + style = BIN or MULTI +------------------------------------------------------------------------- */ + +void NeighList::stencil_allocate(int smax, int style) +{ + int i; + + if (style == BIN) { + if (smax > maxstencil) { + maxstencil = smax; + memory->sfree(stencil); + stencil = (int *) memory->smalloc(smax*sizeof(int),"neighlist:stencil"); + } + + } else { + int n = atom->ntypes; + if (maxstencil_multi == 0) { + nstencil_multi = new int[n+1]; + stencil_multi = new int*[n+1]; + distsq_multi = new double*[n+1]; + for (i = 1; i <= n; i++) { + nstencil_multi[i] = 0; + stencil_multi[i] = NULL; + distsq_multi[i] = NULL; + } + } + if (smax > maxstencil_multi) { + maxstencil_multi = smax; + for (i = 1; i <= n; i++) { + memory->sfree(stencil_multi[i]); + memory->sfree(distsq_multi[i]); + stencil_multi[i] = (int *) + memory->smalloc(smax*sizeof(int),"neighlist:stencil_multi"); + distsq_multi[i] = (double *) memory->smalloc(smax*sizeof(double), + "neighlist:distsq_multi"); + } + } + } +} + +/* ---------------------------------------------------------------------- + add PGDELTA pages to neighbor list +------------------------------------------------------------------------- */ + +int **NeighList::add_pages() +{ + int npage = maxpage; + maxpage += PGDELTA; + + pages = (int **) + memory->srealloc(pages,maxpage*sizeof(int *),"neighlist:pages"); + for (int i = npage; i < maxpage; i++) + pages[i] = (int *) memory->smalloc(pgsize*sizeof(int), + "neighlist:pages[i]"); + + if (dnum) { + dpages = (double **) + memory->srealloc(dpages,maxpage*sizeof(double *),"neighlist:dpages"); + for (int i = npage; i < maxpage; i++) + dpages[i] = (double *) memory->smalloc(dnum*pgsize*sizeof(double), + "neighlist:dpages[i]"); + } + + return pages; +} + +/* ---------------------------------------------------------------------- + print attributes of this list and associated request +------------------------------------------------------------------------- */ + +void NeighList::print_attributes() +{ + if (comm->me != 0) return; + + NeighRequest *rq = neighbor->requests[index]; + + printf("Neighbor list/request %d:\n",index); + printf(" %d = build flag\n",buildflag); + printf(" %d = grow flag\n",growflag); + printf(" %d = stencil flag\n",stencilflag); + printf("\n"); + printf(" %d = pair\n",rq->pair); + printf(" %d = fix\n",rq->fix); + printf(" %d = compute\n",rq->compute); + printf(" %d = command\n",rq->command); + printf("\n"); + printf(" %d = half\n",rq->half); + printf(" %d = full\n",rq->full); + printf(" %d = gran\n",rq->gran); + printf(" %d = granhistory\n",rq->granhistory); + printf(" %d = respainner\n",rq->respainner); + printf(" %d = respamiddle\n",rq->respamiddle); + printf(" %d = respaouter\n",rq->respaouter); + printf(" %d = half_from_full\n",rq->half_from_full); + printf("\n"); + printf(" %d = occasional\n",rq->occasional); + printf(" %d = dnum\n",rq->dnum); + printf(" %d = copy\n",rq->copy); + printf(" %d = skip\n",rq->skip); + printf(" %d = otherlist\n",rq->otherlist); + printf(" %p = listskip\n",listskip); + printf("\n"); +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory + if growflag = 0, maxlocal & maxpage will also be 0 + if stencilflag = 0, maxstencil * maxstencil_multi will also be 0 +------------------------------------------------------------------------- */ + +int NeighList::memory_usage() +{ + int bytes = 0; + + bytes += 2 * maxlocal * sizeof(int); + bytes += maxlocal * sizeof(int *); + bytes += maxpage*pgsize * sizeof(int); + + if (dnum) { + bytes += maxlocal * sizeof(double *); + bytes += dnum * maxpage*pgsize * sizeof(double); + } + + if (maxstencil) bytes += maxstencil * sizeof(int); + if (maxstencil_multi) { + bytes += atom->ntypes * maxstencil_multi * sizeof(int); + bytes += atom->ntypes * maxstencil_multi * sizeof(double); + } + + return bytes; +} diff --git a/src/neigh_list.h b/src/neigh_list.h new file mode 100644 index 0000000000..188c701b82 --- /dev/null +++ b/src/neigh_list.h @@ -0,0 +1,88 @@ +/* ---------------------------------------------------------------------- + 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 NEIGH_LIST_H +#define NEIGH_LIST_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +class NeighList : protected Pointers { + public: + int index; // index of which neigh list it is + // needed when a class invokes it directly + + int buildflag; // 1 if pair_build invoked every reneigh + int growflag; // 1 if stores atom-based arrays & pages + int stencilflag; // 1 if stores stencil arrays + + // data structs to store neighbor pairs I,J and associated values + + int inum; // # of I atoms neighbors are stored for + int *ilist; // local indices of I atoms + int *numneigh; // # of J neighbors for each I atom + int **firstneigh; // ptr to 1st J int value of each I atom + double **firstdouble; // ptr to 1st J double value of each I atom + + int pgsize; // size of each page + int maxpage; // # of pages currently allocated + int **pages; // neighbor list pages for ints + double **dpages; // neighbor list pages for doubles + int dnum; // # of doubles for each pair (0 if none) + + // atom types to skip when building list + // iskip,ijskip are just ptrs to corresponding request + + int skip; // 1 if this list skips atom types from another list + int *iskip; // iskip[i] if atoms of type I are not in list + int **ijskip; // ijskip[i][j] if pairs of type I,J are not in list + + // settings and pointers for related neighbor lists and fixes + + NeighList *listgranhistory; // point at history list + class FixShearHistory *fix_history; // fix that stores history info + + int respamiddle; // 1 if this respaouter has middle list + NeighList *listinner; // me = respaouter, point to respainner + NeighList *listmiddle; // me = respaouter, point to respamiddle + NeighList *listfull; // me = half list, point to full I derive from + NeighList *listcopy; // me = copy list, point to list I copy from + NeighList *listskip; // me = skip list, point to list I skip from + + // stencils of bin indices for neighbor finding + + int maxstencil; // max size of stencil + int nstencil; // # of bins in stencil + int *stencil; // list of bin offsets + + int maxstencil_multi; // max sizes of stencils + int *nstencil_multi; // # bins in each type-based multi stencil + int **stencil_multi; // list of bin offsets in each stencil + double **distsq_multi; // sq distances to bins in each stencil + + NeighList(class LAMMPS *, int); + ~NeighList(); + void grow(int); // grow maxlocal + void stencil_allocate(int, int); // allocate stencil arrays + int **add_pages(); // add pages to neigh list + void print_attributes(); // debug routine + int memory_usage(); + + private: + int maxlocal; // size of allocated atom arrays +}; + +} + +#endif diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp new file mode 100644 index 0000000000..5739a8175c --- /dev/null +++ b/src/neigh_request.cpp @@ -0,0 +1,164 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "neigh_request.h" +#include "atom.h" +#include "memory.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp) +{ + // default ID = 0 + + id = 0; + + // default is pair request + + pair = 1; + fix = compute = command = 0; + + // default is half neighbor list + + half = 1; + full = 0; + gran = granhistory = 0; + respainner = respamiddle = respaouter = 0; + half_from_full = 0; + + // default is every reneighboring + + occasional = 0; + + // default is no auxiliary floating point values + + dnum = 0; + + // default is no copy or skip + + copy = 0; + skip = 0; + iskip = NULL; + ijskip = NULL; + otherlist = -1; +} + +/* ---------------------------------------------------------------------- */ + +NeighRequest::~NeighRequest() +{ + delete [] iskip; + memory->destroy_2d_int_array(ijskip); +} + +/* ---------------------------------------------------------------------- + compare this request to other request + return 1 if every variable is same, 0 if different in any way +------------------------------------------------------------------------- */ + +int NeighRequest::identical(NeighRequest *other) +{ + int same = 1; + + if (requestor != other->requestor) same = 0; + if (id != other->id) same = 0; + + if (pair != other->pair) same = 0; + if (fix != other->fix) same = 0; + if (compute != other->compute) same = 0; + if (command != other->command) same = 0; + + if (half != other->half) same = 0; + if (full != other->full) same = 0; + if (gran != other->gran) same = 0; + if (granhistory != other->granhistory) same = 0; + if (respainner != other->respainner) same = 0; + if (respamiddle != other->respamiddle) same = 0; + if (respaouter != other->respaouter) same = 0; + if (half_from_full != other->half_from_full) same = 0; + + if (occasional != other->occasional) same = 0; + if (dnum != other->dnum) same = 0; + + if (copy != other->copy) same = 0; + if (same_skip(other) == 0) same = 0; + if (otherlist != other->otherlist) same = 0; + + return same; +} + +/* ---------------------------------------------------------------------- + compare kind of this request to other request + return 1 if same, 0 if different +------------------------------------------------------------------------- */ + +int NeighRequest::same_kind(NeighRequest *other) +{ + int same = 1; + + if (half != other->half) same = 0; + if (full != other->full) same = 0; + if (gran != other->gran) same = 0; + if (granhistory != other->granhistory) same = 0; + if (respainner != other->respainner) same = 0; + if (respamiddle != other->respamiddle) same = 0; + if (respaouter != other->respaouter) same = 0; + if (half_from_full != other->half_from_full) same = 0; + + return same; +} + +/* ---------------------------------------------------------------------- + set kind value of this request to that of other request +------------------------------------------------------------------------- */ + +void NeighRequest::copy_kind(NeighRequest *other) +{ + half = 0; + + if (other->half) half = 1; + if (other->full) full = 1; + if (other->gran) gran = 1; + if (other->granhistory) granhistory = 1; + if (other->respainner) respainner = 1; + if (other->respamiddle) respamiddle = 1; + if (other->respaouter) respaouter = 1; + if (other->half_from_full) half_from_full = 1; +} + +/* ---------------------------------------------------------------------- + compare skip attributes of this request to other request + return 1 if same, 0 if different +------------------------------------------------------------------------- */ + +int NeighRequest::same_skip(NeighRequest *other) +{ + int i,j; + + int same = 1; + + if (skip != other->skip) same = 0; + if (skip && other->skip) { + int ntypes = atom->ntypes; + for (i = 1; i <= ntypes; i++) + if (iskip[i] != other->iskip[i]) same = 0; + for (i = 1; i <= ntypes; i++) + for (j = 1; j <= ntypes; j++) + if (ijskip[i][j] != other->ijskip[i][j]) same = 0; + } + + return same; +} + diff --git a/src/neigh_request.h b/src/neigh_request.h new file mode 100644 index 0000000000..9c0a43af60 --- /dev/null +++ b/src/neigh_request.h @@ -0,0 +1,83 @@ +/* ---------------------------------------------------------------------- + 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 NEIGH_REQUEST_H +#define NEIGH_REQUEST_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +class NeighRequest : protected Pointers { + public: + void *requestor; // class that made request + int id; // ID of request + // used to track multiple requests from one class + + // which class is requesting the list, one flag is 1, others are 0 + + int pair; + int fix; + int compute; + int command; + + // kind of list requested, one flag is 1, others are 0 + // set by reqeusting class + + int half; // 1 if half neigh list + int full; // 1 if full neigh list + + int gran; // 1 if granular list + int granhistory; // 1 if granular history list + + int respainner; // 1 if a rRESPA inner list + int respamiddle; // 1 if a rRESPA middle list + int respaouter; // 1 if a rRESPA outer list + + int half_from_full; // 1 if half list computed from previous full list + + // 0 if needed every reneighboring during run + // 1 if occasionally needed by a fix, compute, etc + // set by requesting class + + int occasional; + + // number of auxiliary floating point values to store, 0 if none + // set by requesting class + + int dnum; + + // set by neighbor and pair_hybrid after all requests are made + // these settings do not change kind value + + int copy; // 1 if this list copied from another list + + int skip; // 1 if this list skips atom types from another list + int *iskip; // iskip[i] if atoms of type I are not in list + int **ijskip; // ijskip[i][j] if pairs of type I,J are not in list + + int otherlist; // other list to copy or skip from + + // methods + + NeighRequest(class LAMMPS *); + ~NeighRequest(); + int identical(NeighRequest *); + int same_kind(NeighRequest *); + void copy_kind(NeighRequest *); + int same_skip(NeighRequest *); +}; + +} + +#endif diff --git a/src/pair_hybrid_overlay.cpp b/src/pair_hybrid_overlay.cpp new file mode 100644 index 0000000000..501018aa2f --- /dev/null +++ b/src/pair_hybrid_overlay.cpp @@ -0,0 +1,140 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "pair_hybrid_overlay.h" +#include "atom.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_request.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairHybridOverlay::PairHybridOverlay(LAMMPS *lmp) : PairHybrid(lmp) {} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairHybridOverlay::coeff(int narg, char **arg) +{ + if (narg < 3) error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + // 3rd arg = pair sub-style name + // allow for "none" as valid sub-style name + + int m; + for (m = 0; m < nstyles; m++) + if (strcmp(arg[2],keywords[m]) == 0) break; + + int none = 0; + if (m == nstyles) { + if (strcmp(arg[2],"none") == 0) none = 1; + else error->all("Pair coeff for hybrid has invalid style"); + } + + // move 1st/2nd args to 2nd/3rd args + + sprintf(arg[2],"%s",arg[1]); + sprintf(arg[1],"%s",arg[0]); + + // invoke sub-style coeff() starting with 1st arg + + if (!none) styles[m]->coeff(narg-1,&arg[1]); + + // set setflag and which type pairs map to which sub-style + // if sub-style is none: set hybrid subflag, wipe out map + // else: set hybrid setflag & map only if substyle setflag is set + // multiple mappings are allowed, + // but don't map again if already mapped to this sub-style + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + if (none) { + setflag[i][j] = 1; + nmap[i][j] = 0; + count++; + } else if (styles[m]->setflag[i][j] ) { + int k; + for (k = 0; k < nmap[i][j]; k++) + if (map[i][j][k] == m) break; + if (k < nmap[i][j]) continue; + setflag[i][j] = 1; + map[i][j][nmap[i][j]++] = m; + count++; + } + } + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + combine sub-style neigh list requests and create new ones if needed +------------------------------------------------------------------------- */ + +void PairHybridOverlay::modify_requests() +{ + int i,j; + NeighRequest *irq,*jrq; + + // loop over pair requests + // if a previous list is same kind with same skip attributes + // then make this one a copy list of that one + // NOTE: SHOULD I SKIP granhistory, respamiddle, halffromfull ?? + + for (i = 0; i < neighbor->nrequest; i++) { + irq = neighbor->requests[i]; + for (j = 0; j < i; j++) { + jrq = neighbor->requests[j]; + if (irq->same_kind(jrq) && irq->same_skip(jrq)) { + irq->copy = 1; + irq->otherlist = j; + break; + } + } + } + + // if list is skip list and not copy, look for non-skip list of same kind + // if one exists, point at that one + // else make new non-skip request of same kind and point at that one + + for (i = 0; i < neighbor->nrequest; i++) { + irq = neighbor->requests[i]; + if (irq->skip == 0 || irq->copy) continue; + + for (j = 0; j < neighbor->nrequest; j++) { + jrq = neighbor->requests[j]; + if (irq->same_kind(jrq) && jrq->skip == 0) break; + } + + if (j < neighbor->nrequest) irq->otherlist = j; + else { + int newrequest = neighbor->request(this); + neighbor->requests[newrequest]->copy_kind(irq); + irq->otherlist = newrequest; + } + } +} diff --git a/src/pair_hybrid_overlay.h b/src/pair_hybrid_overlay.h new file mode 100644 index 0000000000..4603f47bc5 --- /dev/null +++ b/src/pair_hybrid_overlay.h @@ -0,0 +1,32 @@ +/* ---------------------------------------------------------------------- + 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 PAIR_HYBRID_OVERLAY_H +#define PAIR_HYBRID_OVERLAY_H + +#include "pair_hybrid.h" + +namespace LAMMPS_NS { + +class PairHybridOverlay : public PairHybrid { + public: + PairHybridOverlay(class LAMMPS *); + void coeff(int, char **); + + private: + void modify_requests(); +}; + +} + +#endif