diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp new file mode 100644 index 0000000000..b57cdc049b --- /dev/null +++ b/src/compute_angle_local.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 "math.h" +#include "string.h" +#include "compute_angle_local.h" +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "domain.h" +#include "force.h" +#include "angle.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{THETA,ENERGY}; + +#define DELTA 10000 + +/* ---------------------------------------------------------------------- */ + +ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 4) error->all("Illegal compute angle/local command"); + + if (atom->avec->angles_allow == 0) + error->all("Compute angle/local used when angles are not allowed"); + + local_flag = 1; + nvalues = narg - 3; + if (nvalues == 1) size_local_cols = 0; + else size_local_cols = nvalues; + + which = new int[nvalues]; + pack_choice = new FnPtrPack[nvalues]; + + dflag = eflag = 0; + + int i; + for (int iarg = 3; iarg < narg; iarg++) { + i = iarg-3; + + if (strcmp(arg[iarg],"theta") == 0) { + dflag = 1; + which[i] = THETA; + pack_choice[i] = &ComputeAngleLocal::pack_theta; + } else if (strcmp(arg[iarg],"energy") == 0) { + eflag = 1; + which[i] = ENERGY; + pack_choice[i] = &ComputeAngleLocal::pack_energy; + } else error->all("Invalid keyword in compute angle/local command"); + } + + nmax = 0; + theta = energy = NULL; + array = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeAngleLocal::~ComputeAngleLocal() +{ + delete [] which; + delete [] pack_choice; + memory->sfree(theta); + memory->sfree(energy); + memory->destroy_2d_double_array(array); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeAngleLocal::init() +{ + if (force->angle == NULL) + error->all("No angle style is defined for compute angle/local"); + + // do initial memory allocation so that memory_usage() is correct + + ncount = compute_angles(0); + if (ncount > nmax) reallocate(ncount); + size_local_rows = ncount; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeAngleLocal::compute_local() +{ + invoked_local = update->ntimestep; + + // count local entries and compute angle info + + ncount = compute_angles(0); + if (ncount > nmax) reallocate(ncount); + size_local_rows = ncount; + ncount = compute_angles(1); + + // fill array with theta/energy values + + if (nvalues > 1) { + if (array) buf = array[0]; + for (int n = 0; n < nvalues; n++) + (this->*pack_choice[n])(n); + } +} + +/* ---------------------------------------------------------------------- + count angles and compute angle info on this proc + only count angle once if newton_angle is off + all atoms in interaction must be in group + all atoms in interaction must be known to proc + if angle is deleted (type = 0), do not count + if angle is turned off (type < 0), still count + if flag is set, compute requested info about angle +------------------------------------------------------------------------- */ + +int ComputeAngleLocal::compute_angles(int flag) +{ + int i,atom1,atom2,atom3; + double delx1,dely1,delz1,delx2,dely2,delz2; + double rsq1,rsq2,r1,r2,c; + + double **x = atom->x; + int *num_angle = atom->num_angle; + int **angle_atom1 = atom->angle_atom1; + int **angle_atom2 = atom->angle_atom2; + int **angle_atom3 = atom->angle_atom3; + int **angle_type = atom->angle_type; + int *tag = atom->tag; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + Angle *angle = force->angle; + + int m = 0; + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + for (i = 0; i < num_angle[atom2]; i++) { + if (tag[atom2] != angle_atom2[atom2][i]) continue; + atom1 = atom->map(angle_atom1[atom2][i]); + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + atom3 = atom->map(angle_atom3[atom2][i]); + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + + if (flag) { + + if (dflag) { + delx1 = x[atom1][0] - x[atom2][0]; + dely1 = x[atom1][1] - x[atom2][1]; + delz1 = x[atom1][2] - x[atom2][2]; + domain->minimum_image(delx1,dely1,delz1); + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + r1 = sqrt(rsq1); + + delx2 = x[atom3][0] - x[atom2][0]; + dely2 = x[atom3][1] - x[atom2][1]; + delz2 = x[atom3][2] - x[atom2][2]; + domain->minimum_image(delx2,dely2,delz2); + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // c = cosine of angle + + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + theta[m] = acos(c); + } + + if (eflag) + energy[m] = angle->single(angle_type[atom1][i],atom1,atom2,atom3); + } + + m++; + } + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeAngleLocal::pack_theta(int n) +{ + for (int m = 0; m < ncount; m++) { + buf[n] = theta[m]; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeAngleLocal::pack_energy(int n) +{ + for (int m = 0; m < ncount; m++) { + buf[n] = energy[m]; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeAngleLocal::reallocate(int n) +{ + // grow vector or array and indices array + + while (nmax < n) nmax += DELTA; + + if (dflag) { + memory->sfree(theta); + theta = (double *) memory->smalloc(nmax*sizeof(double), + "bond/local:theta"); + } + if (eflag) { + memory->sfree(energy); + energy = (double *) memory->smalloc(nmax*sizeof(double), + "bond/local:energy"); + } + + if (nvalues == 1) { + if (dflag) vector_local = theta; + if (eflag) vector_local = energy; + } else { + memory->destroy_2d_double_array(array); + array = memory->create_2d_double_array(nmax,nvalues, + "bond/local:array"); + array_local = array; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeAngleLocal::memory_usage() +{ + double bytes = 0.0; + if (dflag) bytes += nmax * sizeof(double); + if (eflag) bytes += nmax * sizeof(double); + if (nvalues > 1) bytes += nmax*nvalues * sizeof(double); + return bytes; +} diff --git a/src/compute_angle_local.h b/src/compute_angle_local.h new file mode 100644 index 0000000000..bc1bf849f4 --- /dev/null +++ b/src/compute_angle_local.h @@ -0,0 +1,52 @@ +/* ---------------------------------------------------------------------- + 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_ANGLE_LOCAL_H +#define COMPUTE_ANGLE_LOCAL_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeAngleLocal : public Compute { + public: + ComputeAngleLocal(class LAMMPS *, int, char **); + ~ComputeAngleLocal(); + void init(); + void compute_local(); + double memory_usage(); + + private: + int nvalues,dflag,eflag; + int *which; + int ncount; + + int nmax; + double *theta; + double *energy; + double **array; + double *buf; + + int compute_angles(int); + void reallocate(int); + + typedef void (ComputeAngleLocal::*FnPtrPack)(int); + FnPtrPack *pack_choice; // ptrs to pack functions + + void pack_theta(int); + void pack_energy(int); +}; + +} + +#endif diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp new file mode 100644 index 0000000000..a0179d29c2 --- /dev/null +++ b/src/compute_bond_local.cpp @@ -0,0 +1,233 @@ +/* ---------------------------------------------------------------------- + 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 "math.h" +#include "string.h" +#include "compute_bond_local.h" +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "domain.h" +#include "force.h" +#include "bond.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{DISTANCE,ENERGY}; + +#define DELTA 10000 + +/* ---------------------------------------------------------------------- */ + +ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 4) error->all("Illegal compute bond/local command"); + + if (atom->avec->bonds_allow == 0) + error->all("Compute bond/local used when bonds are not allowed"); + + local_flag = 1; + nvalues = narg - 3; + if (nvalues == 1) size_local_cols = 0; + else size_local_cols = nvalues; + + which = new int[nvalues]; + pack_choice = new FnPtrPack[nvalues]; + + dflag = eflag = 0; + + int i; + for (int iarg = 3; iarg < narg; iarg++) { + i = iarg-3; + + if (strcmp(arg[iarg],"distance") == 0) { + dflag = 1; + which[i] = DISTANCE; + pack_choice[i] = &ComputeBondLocal::pack_distance; + } else if (strcmp(arg[iarg],"energy") == 0) { + eflag = 1; + which[i] = ENERGY; + pack_choice[i] = &ComputeBondLocal::pack_energy; + } else error->all("Invalid keyword in compute bond/local command"); + } + + nmax = 0; + distance = energy = NULL; + array = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeBondLocal::~ComputeBondLocal() +{ + delete [] which; + delete [] pack_choice; + memory->sfree(distance); + memory->sfree(energy); + memory->destroy_2d_double_array(array); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBondLocal::init() +{ + if (force->bond == NULL) + error->all("No bond style is defined for compute bond/local"); + + // do initial memory allocation so that memory_usage() is correct + + ncount = compute_bonds(0); + if (ncount > nmax) reallocate(ncount); + size_local_rows = ncount; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBondLocal::compute_local() +{ + invoked_local = update->ntimestep; + + // count local entries and compute bond info + + ncount = compute_bonds(0); + if (ncount > nmax) reallocate(ncount); + size_local_rows = ncount; + ncount = compute_bonds(1); + + // fill array with distance/energy values + + if (nvalues > 1) { + if (array) buf = array[0]; + for (int n = 0; n < nvalues; n++) + (this->*pack_choice[n])(n); + } +} + +/* ---------------------------------------------------------------------- + count bonds and compute bond info on this proc + only count bond once if newton_bond is off + all atoms in interaction must be in group + all atoms in interaction must be known to proc + if bond is deleted (type = 0), do not count + if bond is turned off (type < 0), still count + if flag is set, compute requested info about bond +------------------------------------------------------------------------- */ + +int ComputeBondLocal::compute_bonds(int flag) +{ + int i,atom1,atom2; + double delx,dely,delz,rsq; + + double **x = atom->x; + int *num_bond = atom->num_bond; + int **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *tag = atom->tag; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + Bond *bond = force->bond; + + int m = 0; + for (atom1 = 0; atom1 < nlocal; atom1++) { + if (!(mask[atom1] & groupbit)) continue; + for (i = 0; i < num_bond[atom1]; i++) { + atom2 = atom->map(bond_atom[atom1][i]); + if (atom2 < 0 || !(mask[atom2] & groupbit)) continue; + if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue; + if (bond_type[atom1][i] == 0) continue; + + if (flag) { + delx = x[atom1][0] - x[atom2][0]; + dely = x[atom1][1] - x[atom2][1]; + delz = x[atom1][2] - x[atom2][2]; + domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + if (dflag) distance[m] = sqrt(rsq); + if (eflag) + energy[m] = bond->single(bond_type[atom1][i],rsq,atom1,atom2); + } + + m++; + } + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBondLocal::pack_distance(int n) +{ + for (int m = 0; m < ncount; m++) { + buf[n] = distance[m]; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBondLocal::pack_energy(int n) +{ + for (int m = 0; m < ncount; m++) { + buf[n] = energy[m]; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBondLocal::reallocate(int n) +{ + // grow vector or array and indices array + + while (nmax < n) nmax += DELTA; + + if (dflag) { + memory->sfree(distance); + distance = (double *) memory->smalloc(nmax*sizeof(double), + "bond/local:distance"); + } + if (eflag) { + memory->sfree(energy); + energy = (double *) memory->smalloc(nmax*sizeof(double), + "bond/local:energy"); + } + + if (nvalues == 1) { + if (dflag) vector_local = distance; + if (eflag) vector_local = energy; + } else { + memory->destroy_2d_double_array(array); + array = memory->create_2d_double_array(nmax,nvalues, + "bond/local:array"); + array_local = array; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeBondLocal::memory_usage() +{ + double bytes = 0.0; + if (dflag) bytes += nmax * sizeof(double); + if (eflag) bytes += nmax * sizeof(double); + if (nvalues > 1) bytes += nmax*nvalues * sizeof(double); + return bytes; +} diff --git a/src/compute_bond_local.h b/src/compute_bond_local.h new file mode 100644 index 0000000000..8b6034e330 --- /dev/null +++ b/src/compute_bond_local.h @@ -0,0 +1,52 @@ +/* ---------------------------------------------------------------------- + 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_BOND_LOCAL_H +#define COMPUTE_BOND_LOCAL_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeBondLocal : public Compute { + public: + ComputeBondLocal(class LAMMPS *, int, char **); + ~ComputeBondLocal(); + void init(); + void compute_local(); + double memory_usage(); + + private: + int nvalues,dflag,eflag; + int *which; + int ncount; + + int nmax; + double *distance; + double *energy; + double **array; + double *buf; + + int compute_bonds(int); + void reallocate(int); + + typedef void (ComputeBondLocal::*FnPtrPack)(int); + FnPtrPack *pack_choice; // ptrs to pack functions + + void pack_distance(int); + void pack_energy(int); +}; + +} + +#endif diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index 808046936c..76daabf825 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -40,7 +40,7 @@ enum{NONE,BOND,ANGLE,DIHEDRAL,IMPROPER}; ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - if (narg < 4) error->all("Illegal compute property/atom command"); + if (narg < 4) error->all("Illegal compute property/local command"); local_flag = 1; nvalues = narg - 3; @@ -233,6 +233,8 @@ void ComputePropertyLocal::compute_local() int ComputePropertyLocal::count_bonds(int flag) { + int i,atom1,atom2; + int *num_bond = atom->num_bond; int **bond_atom = atom->bond_atom; int **bond_type = atom->bond_type; @@ -241,8 +243,6 @@ int ComputePropertyLocal::count_bonds(int flag) int nlocal = atom->nlocal; int newton_bond = force->newton_bond; - int i,atom1,atom2; - int m = 0; for (atom1 = 0; atom1 < nlocal; atom1++) { if (!(mask[atom1] & groupbit)) continue; @@ -272,6 +272,8 @@ int ComputePropertyLocal::count_bonds(int flag) int ComputePropertyLocal::count_angles(int flag) { + int i,atom1,atom2,atom3; + int *num_angle = atom->num_angle; int **angle_atom1 = atom->angle_atom1; int **angle_atom2 = atom->angle_atom2; @@ -281,8 +283,6 @@ int ComputePropertyLocal::count_angles(int flag) int *mask = atom->mask; int nlocal = atom->nlocal; - int i,atom1,atom2,atom3; - int m = 0; for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; @@ -313,6 +313,8 @@ int ComputePropertyLocal::count_angles(int flag) int ComputePropertyLocal::count_dihedrals(int flag) { + int i,atom1,atom2,atom3,atom4; + int *num_dihedral = atom->num_dihedral; int **dihedral_atom1 = atom->dihedral_atom1; int **dihedral_atom2 = atom->dihedral_atom2; @@ -323,8 +325,6 @@ int ComputePropertyLocal::count_dihedrals(int flag) int *mask = atom->mask; int nlocal = atom->nlocal; - int i,atom1,atom2,atom3,atom4; - int m = 0; for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; @@ -357,6 +357,8 @@ int ComputePropertyLocal::count_dihedrals(int flag) int ComputePropertyLocal::count_impropers(int flag) { + int i,atom1,atom2,atom3,atom4; + int *num_improper = atom->num_improper; int **improper_atom1 = atom->improper_atom1; int **improper_atom2 = atom->improper_atom2; @@ -367,8 +369,6 @@ int ComputePropertyLocal::count_impropers(int flag) int *mask = atom->mask; int nlocal = atom->nlocal; - int i,atom1,atom2,atom3,atom4; - int m = 0; for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; diff --git a/src/style.h b/src/style.h index ce4aa7da52..ead26a2cb2 100644 --- a/src/style.h +++ b/src/style.h @@ -76,6 +76,8 @@ CommandStyle(write_restart,WriteRestart) #endif #ifdef ComputeInclude +#include "compute_angle_local.h" +#include "compute_bond_local.h" #include "compute_centro_atom.h" #include "compute_cna_atom.h" #include "compute_com.h" @@ -108,6 +110,8 @@ CommandStyle(write_restart,WriteRestart) #endif #ifdef ComputeClass +ComputeStyle(angle/local,ComputeAngleLocal) +ComputeStyle(bond/local,ComputeBondLocal) ComputeStyle(centro/atom,ComputeCentroAtom) ComputeStyle(cna/atom,ComputeCNAAtom) ComputeStyle(com,ComputeCOM)