jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11328 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2014-01-25 22:46:08 +00:00
parent 2f6a46a2be
commit 0ce16af78b
63 changed files with 2943 additions and 798 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/GPU/fix_gpu.cpp
View File

@ -163,6 +163,12 @@ int FixGPU::setmask()
void FixGPU::init()
{
// GPU package cannot be used with atom_style template
if (atom->molecular == 2)
error->all(FLERR,"GPU package does not (yet) work with "
"atom_style template");
// hybrid cannot be used with force/neigh option
if (_gpu_mode == GPU_NEIGH || _gpu_mode == GPU_HYB_NEIGH)
@ -188,7 +194,7 @@ void FixGPU::init()
force->pair->no_virial_fdotr_compute = 1;
}
// r-RESPA support
// rRESPA support
if (strstr(update->integrate_style,"respa"))
_nlevels_respa = ((Respa *) update->integrate)->nlevels;

22
src/GRANULAR/fix_pour.cpp
View File

@ -121,9 +121,14 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Fix pour molecule must have coordinates");
if (onemol->typeflag == 0)
error->all(FLERR,"Fix pour molecule must have atom types");
if (ntype+onemol->maxtype <= 0 || ntype+onemol->maxtype > atom->ntypes)
if (ntype+onemol->ntypes <= 0 || ntype+onemol->ntypes > atom->ntypes)
error->all(FLERR,"Invalid atom type in fix pour mol command");
if (atom->molecular == 2 && onemol != atom->avec->onemols[0])
error->all(FLERR,"Fix pour molecule template ID must be same "
"as atom style template ID");
onemol->check_attributes(0);
// fix pour uses geoemetric center of molecule for insertion
onemol->compute_center();
@ -325,7 +330,8 @@ void FixPour::init()
int tmp;
if (onemol != (Molecule *) fixrigid->extract("onemol",tmp))
error->all(FLERR,
"Fix pour and fix rigid/small not using same molecule ID");
"Fix pour and fix rigid/small not using "
"same molecule template ID");
}
// if shakeflag defined, check for SHAKE fix
@ -338,7 +344,8 @@ void FixPour::init()
fixshake = modify->fix[ifix];
int tmp;
if (onemol != (Molecule *) fixshake->extract("onemol",tmp))
error->all(FLERR,"Fix pour and fix shake not using same molecule ID");
error->all(FLERR,"Fix pour and fix shake not using "
"same molecule template ID");
}
}
@ -585,8 +592,13 @@ void FixPour::pre_exchange()
else atom->avec->create_atom(ntype+onemol->type[m],coords[m]);
int n = atom->nlocal - 1;
atom->tag[n] = maxtag_all + m+1;
if (mode == MOLECULE && atom->molecule_flag)
atom->molecule[n] = maxmol_all+1;
if (mode == MOLECULE) {
if (atom->molecular) atom->molecule[n] = maxmol_all+1;
if (atom->molecular == 2) {
atom->molindex[n] = 0;
atom->molatom[n] = m;
}
}
atom->mask[n] = 1 | groupbit;
atom->image[n] = imageflags[m];
atom->v[n][0] = vnew[0];

2
src/MC/fix_bond_break.cpp
View File

@ -78,7 +78,7 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) :
// error check
if (atom->molecular == 0)
if (atom->molecular != 1)
error->all(FLERR,"Cannot use fix bond/break with non-molecular systems");
// initialize Marsaglia RNG with processor-unique seed

2
src/MC/fix_bond_create.cpp
View File

@ -107,7 +107,7 @@ FixBondCreate::FixBondCreate(LAMMPS *lmp, int narg, char **arg) :
// error check
if (atom->molecular == 0)
if (atom->molecular != 1)
error->all(FLERR,"Cannot use fix bond/create with non-molecular systems");
if (iatomtype == jatomtype &&
((imaxbond != jmaxbond) || (inewtype != jnewtype)))

5
src/MC/fix_bond_swap.cpp
View File

@ -72,6 +72,11 @@ FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) :
int seed = force->inumeric(FLERR,arg[5]);
random = new RanMars(lmp,seed + comm->me);
// error check
if (atom->molecular != 1)
error->all(FLERR,"Cannot use fix bond/swap with non-molecular systems");
// create a new compute temp style
// id = fix-ID + temp, compute group = fix group

20
src/MC/fix_gcmc.cpp
View File

@ -49,6 +49,9 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 11) error->all(FLERR,"Illegal fix gcmc command");
if (atom->molecular == 2)
error->all(FLERR,"Fix gcmc does not (yet) work with atom_style template");
vector_flag = 1;
size_vector = 8;
global_freq = 1;
@ -828,7 +831,7 @@ void FixGCMC::attempt_molecule_insertion()
double **x = atom->x;
double **v = atom->v;
imageint *image = atom->image;
int *molecule = atom->molecule;
tagint *molecule = atom->molecule;
tagint *tag = atom->tag;
for (int i = 0; i < natoms_per_molecule; i++) {
k += atom->avec->unpack_exchange(&model_atom_buf[k]);
@ -950,7 +953,7 @@ int FixGCMC::pick_random_gas_atom()
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
int FixGCMC::pick_random_gas_molecule()
tagint FixGCMC::pick_random_gas_molecule()
{
int iwhichglobal = static_cast<int> (ngas*random_equal->uniform());
tagint gas_molecule_id = 0;
@ -1076,21 +1079,12 @@ void FixGCMC::get_model_molecule()
// old_atom = original atom class
// atom = new model atom class
// if old_atom style was hybrid, pass sub-style names to create_avec
Atom *old_atom = atom;
atom = new Atom(lmp);
atom->settings(old_atom);
int nstyles = 0;
char **keywords = NULL;
if (strcmp(old_atom->atom_style,"hybrid") == 0) {
AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) old_atom->avec;
nstyles = avec_hybrid->nstyles;
keywords = avec_hybrid->keywords;
}
atom->create_avec(old_atom->atom_style,nstyles,keywords);
atom->create_avec(old_atom->atom_style,
old_atom->avec->nargcopy,old_atom->avec->argcopy);
// assign atom and topology counts in model atom class from old_atom

4
src/MC/fix_gcmc.h
View File

@ -39,10 +39,10 @@ class FixGCMC : public Fix {
void attempt_molecule_rotation();
void attempt_molecule_deletion();
void attempt_molecule_insertion();
double energy(int, int, int, double *);
double energy(int, int, tagint, double *);
int pick_random_gas_atom();
tagint pick_random_gas_molecule();
double molecule_energy(int);
double molecule_energy(tagint);
void get_rotation_matrix(double, double *);
void get_model_molecule();
void update_gas_atoms_list();

22
src/MISC/fix_deposit.cpp
View File

@ -103,9 +103,14 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Fix deposit molecule must have coordinates");
if (onemol->typeflag == 0)
error->all(FLERR,"Fix deposit molecule must have atom types");
if (ntype+onemol->maxtype <= 0 || ntype+onemol->maxtype > atom->ntypes)
if (ntype+onemol->ntypes <= 0 || ntype+onemol->ntypes > atom->ntypes)
error->all(FLERR,"Invalid atom type in fix deposit mol command");
if (atom->molecular == 2 && onemol != atom->avec->onemols[0])
error->all(FLERR,"Fix deposit molecule template ID must be same "
"as atom style template ID");
onemol->check_attributes(0);
// fix deposit uses geoemetric center of molecule for insertion
onemol->compute_center();
@ -216,7 +221,8 @@ void FixDeposit::init()
int tmp;
if (onemol != (Molecule *) fixrigid->extract("onemol",tmp))
error->all(FLERR,
"Fix deposit and fix rigid/small not using same molecule ID");
"Fix deposit and fix rigid/small not using "
"same molecule template ID");
}
// if shakeflag defined, check for SHAKE fix
@ -229,7 +235,8 @@ void FixDeposit::init()
fixshake = modify->fix[ifix];
int tmp;
if (onemol != (Molecule *) fixshake->extract("onemol",tmp))
error->all(FLERR,"Fix deposit and fix shake not using same molecule ID");
error->all(FLERR,"Fix deposit and fix shake not using "
"same molecule template ID");
}
}
@ -442,8 +449,13 @@ void FixDeposit::pre_exchange()
else atom->avec->create_atom(ntype+onemol->type[m],coords[m]);
n = atom->nlocal - 1;
atom->tag[n] = maxtag_all + m+1;
if (mode == MOLECULE && atom->molecule_flag)
atom->molecule[n] = maxmol_all+1;
if (mode == MOLECULE) {
if (atom->molecular) atom->molecule[n] = maxmol_all+1;
if (atom->molecular == 2) {
atom->molindex[n] = 0;
atom->molatom[n] = m;
}
}
atom->mask[n] = 1 | groupbit;
atom->image[n] = imageflags[m];
atom->v[n][0] = vnew[0];

30
src/MISC/fix_evaporate.cpp
View File

@ -17,6 +17,7 @@
#include "fix_evaporate.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h"
#include "region.h"
@ -234,9 +235,13 @@ void FixEvaporate::pre_exchange()
// keep ndel,ndeltopo,ncount,nall,nbefore current after each mol deletion
} else {
int me,proc,iatom,ndelone,ndelall;
tagint imol;
int me,proc,iatom,ndelone,ndelall,index;
tagint imolecule;
tagint *molecule = atom->molecule;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int molecular = atom->molecular;
Molecule **onemols = atom->avec->onemols;
ndeltopo[0] = ndeltopo[1] = ndeltopo[2] = ndeltopo[3] = 0;
@ -248,23 +253,26 @@ void FixEvaporate::pre_exchange()
if (iwhichglobal >= nbefore && iwhichglobal < nbefore + ncount) {
iwhichlocal = iwhichglobal - nbefore;
iatom = list[iwhichlocal];
imol = molecule[iatom];
imolecule = molecule[iatom];
me = comm->me;
} else me = -1;
// bcast mol ID to delete all atoms from
// if mol ID > 0, delete any atom in molecule and decrement counters
// if mol ID == 0, delete single iatom
// be careful to delete correct # of bond,angle,etc for newton on or off
// logic with ndeltopo is to count # of deleted bonds,angles,etc
// for atom->molecular = 1, do this for each deleted atom in molecule
// for atom->molecular = 2, use Molecule counts for just 1st atom in mol
MPI_Allreduce(&me,&proc,1,MPI_INT,MPI_MAX,world);
MPI_Bcast(&imol,1,MPI_LMP_TAGINT,proc,world);
MPI_Bcast(&imolecule,1,MPI_LMP_TAGINT,proc,world);
ndelone = 0;
for (i = 0; i < nlocal; i++) {
if (imol && molecule[i] == imol) {
if (imolecule && molecule[i] == imolecule) {
mark[i] = 1;
ndelone++;
if (molecular == 1) {
if (atom->avec->bonds_allow) {
if (force->newton_bond) ndeltopo[0] += atom->num_bond[i];
else {
@ -301,6 +309,16 @@ void FixEvaporate::pre_exchange()
}
}
} else {
if (molatom[i] == 0) {
index = molindex[i];
ndeltopo[0] += onemols[index]->nbonds;
ndeltopo[1] += onemols[index]->nangles;
ndeltopo[2] += onemols[index]->ndihedrals;
ndeltopo[3] += onemols[index]->nimpropers;
}
}
} else if (me == proc && i == iatom) {
mark[i] = 1;
ndelone++;

867
src/MOLECULE/atom_vec_template.cpp Normal file
View File

@ -0,0 +1,867 @@
/* ----------------------------------------------------------------------
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 "stdlib.h"
#include "atom_vec_template.h"
#include "atom.h"
#include "molecule.h"
#include "comm.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define DELTA 10000
/* ---------------------------------------------------------------------- */
AtomVecTemplate::AtomVecTemplate(LAMMPS *lmp) : AtomVec(lmp)
{
molecular = 2;
mass_type = 1;
comm_x_only = comm_f_only = 1;
size_forward = 3;
size_reverse = 3;
size_border = 9;
size_velocity = 3;
size_data_atom = 8;
size_data_vel = 4;
xcol_data = 6;
atom->molecule_flag = 1;
}
/* ----------------------------------------------------------------------
process additional arg = molecule template ID
------------------------------------------------------------------------- */
void AtomVecTemplate::process_args(int narg, char **arg)
{
if (narg != 1) error->all(FLERR,"Invalid atom_style template command");
int imol = atom->find_molecule(arg[0]);
if (imol == -1) error->all(FLERR,"Molecule template ID for "
"atom_style template does not exist");
onemols = &atom->molecules[imol];
nset = atom->molecules[imol]->nset;
// set bonds_allow,angles_allow,etc based on the molecules in template set
// similar to how atom_style bond,angle,full set it
for (int i = 0; i < nset; i++) {
if (onemols[i]->bondflag) bonds_allow = 1;
if (onemols[i]->angleflag) angles_allow = 1;
if (onemols[i]->dihedralflag) dihedrals_allow = 1;
if (onemols[i]->improperflag) impropers_allow = 1;
}
// set nbondtypes,nangletypes,etc based on the molecules in template set
// do this here b/c data file will typically not contain these settings
for (int i = 0; i < nset; i++) {
atom->nbondtypes = MAX(atom->nbondtypes,onemols[i]->nbondtypes);
atom->nangletypes = MAX(atom->nangletypes,onemols[i]->nangletypes);
atom->ndihedraltypes = MAX(atom->ndihedraltypes,onemols[i]->ndihedraltypes);
atom->nimpropertypes = MAX(atom->nimpropertypes,onemols[i]->nimpropertypes);
}
}
/* ----------------------------------------------------------------------
grow atom arrays
n = 0 grows arrays by DELTA
n > 0 allocates arrays to size n
------------------------------------------------------------------------- */
void AtomVecTemplate::grow(int n)
{
if (n == 0) nmax += DELTA;
else nmax = n;
atom->nmax = nmax;
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");
type = memory->grow(atom->type,nmax,"atom:type");
mask = memory->grow(atom->mask,nmax,"atom:mask");
image = memory->grow(atom->image,nmax,"atom:image");
x = memory->grow(atom->x,nmax,3,"atom:x");
v = memory->grow(atom->v,nmax,3,"atom:v");
f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
molecule = memory->grow(atom->molecule,nmax,"atom:molecule");
molindex = memory->grow(atom->molindex,nmax,"atom:molindex");
molatom = memory->grow(atom->molatom,nmax,"atom:molatom");
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
}
/* ----------------------------------------------------------------------
reset local array ptrs
------------------------------------------------------------------------- */
void AtomVecTemplate::grow_reset()
{
tag = atom->tag; type = atom->type;
mask = atom->mask; image = atom->image;
x = atom->x; v = atom->v; f = atom->f;
molecule = atom->molecule;
molindex = atom->molindex; molatom = atom->molatom;
}
/* ----------------------------------------------------------------------
copy atom I info to atom J
------------------------------------------------------------------------- */
void AtomVecTemplate::copy(int i, int j, int delflag)
{
tag[j] = tag[i];
type[j] = type[i];
mask[j] = mask[i];
image[j] = image[i];
x[j][0] = x[i][0];
x[j][1] = x[i][1];
x[j][2] = x[i][2];
v[j][0] = v[i][0];
v[j][1] = v[i][1];
v[j][2] = v[i][2];
molecule[j] = molecule[i];
molindex[j] = molindex[i];
molatom[j] = molatom[i];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j,delflag);
}
/* ---------------------------------------------------------------------- */
int AtomVecTemplate::pack_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
dz = pbc[2]*domain->zprd;
}
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
}
}
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecTemplate::pack_comm_vel(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz,dvx,dvy,dvz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
dz = pbc[2]*domain->zprd;
}
if (!deform_vremap) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
dvz = pbc[2]*h_rate[2];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
buf[m++] = v[j][2] + dvz;
} else {
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
}
}
}
return m;
}
/* ---------------------------------------------------------------------- */
void AtomVecTemplate::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
void AtomVecTemplate::unpack_comm_vel(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int AtomVecTemplate::pack_reverse(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = f[i][0];
buf[m++] = f[i][1];
buf[m++] = f[i][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void AtomVecTemplate::unpack_reverse(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
f[j][0] += buf[m++];
f[j][1] += buf[m++];
f[j][2] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int AtomVecTemplate::pack_border(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d;
buf[m++] = ubuf(molindex[j]).d;
buf[m++] = ubuf(molatom[j]).d;
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0];
dy = pbc[1];
dz = pbc[2];
}
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d;
buf[m++] = ubuf(molindex[j]).d;
buf[m++] = ubuf(molatom[j]).d;
}
}
if (atom->nextra_border)
for (int iextra = 0; iextra < atom->nextra_border; iextra++)
m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecTemplate::pack_border_vel(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz,dvx,dvy,dvz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d;
buf[m++] = ubuf(molindex[j]).d;
buf[m++] = ubuf(molatom[j]).d;
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0];
dy = pbc[1];
dz = pbc[2];
}
if (!deform_vremap) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d;
buf[m++] = ubuf(molindex[j]).d;
buf[m++] = ubuf(molatom[j]).d;
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
dvz = pbc[2]*h_rate[2];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = ubuf(molecule[j]).d;
buf[m++] = ubuf(molindex[j]).d;
buf[m++] = ubuf(molatom[j]).d;
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
buf[m++] = v[j][2] + dvz;
} else {
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
}
}
}
if (atom->nextra_border)
for (int iextra = 0; iextra < atom->nextra_border; iextra++)
m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecTemplate::pack_border_hybrid(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = ubuf(molecule[j]).d;
buf[m++] = ubuf(molindex[j]).d;
buf[m++] = ubuf(molatom[j]).d;
}
return m;
}
/* ---------------------------------------------------------------------- */
void AtomVecTemplate::unpack_border(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (i == nmax) grow(0);
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
tag[i] = (tagint) ubuf(buf[m++]).i;
type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i;
molecule[i] = (tagint) ubuf(buf[m++]).i;
molindex[i] = (int) ubuf(buf[m++]).i;
molatom[i] = (int) ubuf(buf[m++]).i;
}
if (atom->nextra_border)
for (int iextra = 0; iextra < atom->nextra_border; iextra++)
m += modify->fix[atom->extra_border[iextra]]->
unpack_border(n,first,&buf[m]);
}
/* ---------------------------------------------------------------------- */
void AtomVecTemplate::unpack_border_vel(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (i == nmax) grow(0);
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
tag[i] = (tagint) ubuf(buf[m++]).i;
type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i;
molecule[i] = (tagint) ubuf(buf[m++]).i;
molindex[i] = (int) ubuf(buf[m++]).i;
molatom[i] = (int) ubuf(buf[m++]).i;
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
}
if (atom->nextra_border)
for (int iextra = 0; iextra < atom->nextra_border; iextra++)
m += modify->fix[atom->extra_border[iextra]]->
unpack_border(n,first,&buf[m]);
}
/* ---------------------------------------------------------------------- */
int AtomVecTemplate::unpack_border_hybrid(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
molecule[i] = (tagint) ubuf(buf[m++]).i;
molindex[i] = (int) ubuf(buf[m++]).i;
molatom[i] = (int) ubuf(buf[m++]).i;
}
return m;
}
/* ----------------------------------------------------------------------
pack data for atom I for sending to another proc
xyz must be 1st 3 values, so comm::exchange() can test on them
------------------------------------------------------------------------- */
int AtomVecTemplate::pack_exchange(int i, double *buf)
{
int m = 1;
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
buf[m++] = v[i][0];
buf[m++] = v[i][1];
buf[m++] = v[i][2];
buf[m++] = ubuf(tag[i]).d;
buf[m++] = ubuf(type[i]).d;
buf[m++] = ubuf(mask[i]).d;
buf[m++] = ubuf(image[i]).d;
buf[m++] = ubuf(molecule[i]).d;
buf[m++] = ubuf(molindex[i]).d;
buf[m++] = ubuf(molatom[i]).d;
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
buf[0] = m;
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecTemplate::unpack_exchange(double *buf)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
int m = 1;
x[nlocal][0] = buf[m++];
x[nlocal][1] = buf[m++];
x[nlocal][2] = buf[m++];
v[nlocal][0] = buf[m++];
v[nlocal][1] = buf[m++];
v[nlocal][2] = buf[m++];
tag[nlocal] = (tagint) ubuf(buf[m++]).i;
type[nlocal] = (int) ubuf(buf[m++]).i;
mask[nlocal] = (int) ubuf(buf[m++]).i;
image[nlocal] = (imageint) ubuf(buf[m++]).i;
molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
molindex[nlocal] = (int) ubuf(buf[m++]).i;
molatom[nlocal] = (int) ubuf(buf[m++]).i;
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
m += modify->fix[atom->extra_grow[iextra]]->
unpack_exchange(nlocal,&buf[m]);
atom->nlocal++;
return m;
}
/* ----------------------------------------------------------------------
size of restart data for all atoms owned by this proc
include extra data stored by fixes
------------------------------------------------------------------------- */
int AtomVecTemplate::size_restart()
{
int i;
int nlocal = atom->nlocal;
int n = 14 * nlocal;
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
for (i = 0; i < nlocal; i++)
n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
return n;
}
/* ----------------------------------------------------------------------
pack atom I's data for restart file including extra quantities
xyz must be 1st 3 values, so that read_restart can test on them
molecular types may be negative, but write as positive
------------------------------------------------------------------------- */
int AtomVecTemplate::pack_restart(int i, double *buf)
{
int m = 1;
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
buf[m++] = ubuf(tag[i]).d;
buf[m++] = ubuf(type[i]).d;
buf[m++] = ubuf(mask[i]).d;
buf[m++] = ubuf(image[i]).d;
buf[m++] = v[i][0];
buf[m++] = v[i][1];
buf[m++] = v[i][2];
buf[m++] = ubuf(molecule[i]).d;
buf[m++] = ubuf(molindex[i]).d;
buf[m++] = ubuf(molatom[i]).d;
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
buf[0] = m;
return m;
}
/* ----------------------------------------------------------------------
unpack data for one atom from restart file including extra quantities
------------------------------------------------------------------------- */
int AtomVecTemplate::unpack_restart(double *buf)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) {
grow(0);
if (atom->nextra_store)
memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra");
}
int m = 1;
x[nlocal][0] = buf[m++];
x[nlocal][1] = buf[m++];
x[nlocal][2] = buf[m++];
tag[nlocal] = (tagint) ubuf(buf[m++]).i;
type[nlocal] = (int) ubuf(buf[m++]).i;
mask[nlocal] = (int) ubuf(buf[m++]).i;
image[nlocal] = (imageint) ubuf(buf[m++]).i;
v[nlocal][0] = buf[m++];
v[nlocal][1] = buf[m++];
v[nlocal][2] = buf[m++];
molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
molindex[nlocal] = (int) ubuf(buf[m++]).i;
molatom[nlocal] = (int) ubuf(buf[m++]).i;
double **extra = atom->extra;
if (atom->nextra_store) {
int size = static_cast<int> (buf[0]) - m;
for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
}
atom->nlocal++;
return m;
}
/* ----------------------------------------------------------------------
create one atom of itype at coord
set other values to defaults
------------------------------------------------------------------------- */
void AtomVecTemplate::create_atom(int itype, double *coord)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
tag[nlocal] = 0;
type[nlocal] = itype;
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
molecule[nlocal] = 0;
molindex[nlocal] = -1;
molatom[nlocal] = -1;
atom->nlocal++;
}
/* ----------------------------------------------------------------------
unpack one line from Atoms section of data file
initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecTemplate::data_atom(double *coord, imageint imagetmp, char **values)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
tag[nlocal] = ATOTAGINT(values[0]);
if (tag[nlocal] <= 0)
error->one(FLERR,"Invalid atom ID in Atoms section of data file");
molecule[nlocal] = ATOTAGINT(values[1]);
molindex[nlocal] = atoi(values[2]) - 1;
molatom[nlocal] = atoi(values[3]) - 1;
if (molindex[nlocal] < 0 || molindex[nlocal] >= nset)
error->one(FLERR,"Invalid template index in Atoms section of data file");
if (molatom[nlocal] < 0 ||
molatom[nlocal] >= onemols[molindex[nlocal]]->natoms)
error->one(FLERR,"Invalid template atom in Atoms section of data file");
type[nlocal] = atoi(values[4]);
if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
error->one(FLERR,"Invalid atom type in Atoms section of data file");
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
image[nlocal] = imagetmp;
mask[nlocal] = 1;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
atom->nlocal++;
}
/* ----------------------------------------------------------------------
unpack hybrid quantities from one line in Atoms section of data file
initialize other atom quantities for this sub-style
------------------------------------------------------------------------- */
int AtomVecTemplate::data_atom_hybrid(int nlocal, char **values)
{
molecule[nlocal] = ATOTAGINT(values[0]);
molindex[nlocal] = atoi(values[1]) - 1;
molatom[nlocal] = atoi(values[2]) - 1;
return 3;
}
/* ----------------------------------------------------------------------
pack atom info for data file including 3 image flags
------------------------------------------------------------------------- */
void AtomVecTemplate::pack_data(double **buf)
{
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
buf[i][0] = ubuf(tag[i]).d;
buf[i][1] = ubuf(molecule[i]).d;
buf[i][2] = ubuf(molindex[i]+1).d;
buf[i][3] = ubuf(molatom[i]+1).d;
buf[i][4] = ubuf(type[i]).d;
buf[i][5] = x[i][0];
buf[i][6] = x[i][1];
buf[i][7] = x[i][2];
buf[i][8] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
buf[i][9] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
buf[i][10] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
}
}
/* ----------------------------------------------------------------------
pack hybrid atom info for data file
------------------------------------------------------------------------- */
int AtomVecTemplate::pack_data_hybrid(int i, double *buf)
{
buf[0] = ubuf(molecule[i]).d;
buf[1] = ubuf(molindex[i]+1).d;
buf[2] = ubuf(molatom[i]+1).d;
return 3;
}
/* ----------------------------------------------------------------------
write atom info to data file including 3 image flags
------------------------------------------------------------------------- */
void AtomVecTemplate::write_data(FILE *fp, int n, double **buf)
{
for (int i = 0; i < n; i++)
fprintf(fp,TAGINT_FORMAT " " TAGINT_FORMAT
" %d %d %d %-1.16e %-1.16e %-1.16e %d %d %d\n",
(tagint) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
(int) ubuf(buf[i][2]).i,(int) ubuf(buf[i][3]).i,
(int) ubuf(buf[i][4]).i,
buf[i][5],buf[i][6],buf[i][7],
(int) ubuf(buf[i][8]).i,(int) ubuf(buf[i][9]).i,
(int) ubuf(buf[i][10]).i);
}
/* ----------------------------------------------------------------------
write hybrid atom info to data file
------------------------------------------------------------------------- */
int AtomVecTemplate::write_data_hybrid(FILE *fp, double *buf)
{
fprintf(fp," " TAGINT_FORMAT " %d %d %d",
(tagint) ubuf(buf[0]).i,(int) ubuf(buf[0]).i,(int) ubuf(buf[0]).i);
return 3;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
------------------------------------------------------------------------- */
bigint AtomVecTemplate::memory_usage()
{
bigint bytes = 0;
if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax);
if (atom->memcheck("type")) bytes += memory->usage(type,nmax);
if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax);
if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
if (atom->memcheck("molindex")) bytes += memory->usage(molindex,nmax);
if (atom->memcheck("molatom")) bytes += memory->usage(molatom,nmax);
return bytes;
}

88
src/MOLECULE/atom_vec_template.h Normal file
View File

@ -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.
------------------------------------------------------------------------- */
#ifdef ATOM_CLASS
AtomStyle(template,AtomVecTemplate)
#else
#ifndef LMP_ATOM_VEC_TEMPLATE_H
#define LMP_ATOM_VEC_TEMPLATE_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecTemplate : public AtomVec {
public:
AtomVecTemplate(class LAMMPS *);
virtual ~AtomVecTemplate() {}
void process_args(int, char **);
void grow(int);
void grow_reset();
void copy(int, int, int);
virtual int pack_comm(int, int *, double *, int, int *);
virtual int pack_comm_vel(int, int *, double *, int, int *);
virtual void unpack_comm(int, int, double *);
virtual void unpack_comm_vel(int, int, double *);
int pack_reverse(int, int, double *);
void unpack_reverse(int, int *, double *);
virtual int pack_border(int, int *, double *, int, int *);
virtual int pack_border_vel(int, int *, double *, int, int *);
int pack_border_hybrid(int, int *, double *);
virtual void unpack_border(int, int, double *);
virtual void unpack_border_vel(int, int, double *);
int unpack_border_hybrid(int, int, double *);
virtual int pack_exchange(int, double *);
virtual int unpack_exchange(double *);
int size_restart();
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double *);
void data_atom(double *, tagint, char **);
int data_atom_hybrid(int, char **);
void pack_data(double **);
int pack_data_hybrid(int, double *);
void write_data(FILE *, int, double **);
int write_data_hybrid(FILE *, double *);
bigint memory_usage();
protected:
int *tag,*type,*mask;
tagint *image;
double **x,**v,**f;
int *molecule,*molindex,*molatom;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Per-processor system is too big
The number of owned atoms plus ghost atoms on a single
processor must fit in 32-bit integer.
E: Invalid atom ID in Atoms section of data file
Atom IDs must be positive integers.
E: Invalid atom type in Atoms section of data file
Atom types must range from 1 to specified # of types.
*/

12
src/MOLECULE/bond_quartic.cpp
View File

@ -233,12 +233,12 @@ void BondQuartic::init_style()
{
if (force->pair == NULL || force->pair->single_enable == 0)
error->all(FLERR,"Pair style does not support bond_style quartic");
if (force->angle)
error->all(FLERR,"Bond style quartic cannot be used with 3,4-body interactions");
if (force->dihedral)
error->all(FLERR,"Bond style quartic cannot be used with 3,4-body interactions");
if (force->improper)
error->all(FLERR,"Bond style quartic cannot be used with 3,4-body interactions");
if (force->angle || force->dihedral || force->improper)
error->all(FLERR,
"Bond style quartic cannot be used with 3,4-body interactions");
if (atom->molecular == 2)
error->all(FLERR,
"Bond style quartic cannot be used with atom style template");
// special bonds must be 1 1 1

48
src/MOLECULE/pair_hbond_dreiding_lj.cpp
View File

@ -21,6 +21,8 @@
#include "string.h"
#include "pair_hbond_dreiding_lj.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
@ -77,7 +79,8 @@ PairHbondDreidingLJ::~PairHbondDreidingLJ()
void PairHbondDreidingLJ::compute(int eflag, int vflag)
{
int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype;
int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype,iatom,imol;
tagint tagprev;
double delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch;
double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2,d;
@ -85,6 +88,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
double r2inv,r10inv;
double switch1,switch2;
int *ilist,*jlist,*numneigh,**firstneigh;
int *tlist;
tagint *klist;
evdwl = ehbond = 0.0;
@ -93,10 +97,15 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
tagint **special = atom->special;
int *type = atom->type;
int **nspecial = atom->nspecial;
int *type = atom->type;
double *special_lj = force->special_lj;
int molecular = atom->molecular;
Molecule **onemols = atom->avec->onemols;
inum = list->inum;
ilist = list->ilist;
@ -113,8 +122,17 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
i = ilist[ii];
itype = type[i];
if (!donor[itype]) continue;
if (molecular == 1) {
klist = special[i];
knum = nspecial[i][0];
} else {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
tlist = onemols[imol]->special[iatom];
knum = onemols[imol]->nspecial[iatom][0];
tagprev = tag[i] - iatom - 1;
}
jlist = firstneigh[i];
jnum = numneigh[i];
@ -132,7 +150,8 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
rsq = delx*delx + dely*dely + delz*delz;
for (kk = 0; kk < knum; kk++) {
k = atom->map(klist[kk]);
if (molecular == 1) k = atom->map(klist[kk]);
else k = atom->map(tlist[kk]+tagprev);
if (k < 0) continue;
ktype = type[k];
m = type2param[itype][jtype][ktype];
@ -454,16 +473,16 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype,
double &fforce)
{
int k,kk,ktype,knum,m;
tagint tagprev;
double eng,eng_lj,force_kernel,force_angle;
double rsq1,rsq2,r1,r2,c,s,ac,r2inv,r10inv,factor_hb;
double switch1,switch2;
double delr1[3],delr2[3];
int *tlist;
tagint *klist;
double **x = atom->x;
tagint **special = atom->special;
int *type = atom->type;
int **nspecial = atom->nspecial;
double *special_lj = force->special_lj;
eng = 0.0;
@ -474,13 +493,26 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype,
if (!donor[itype]) return 0.0;
if (!acceptor[jtype]) return 0.0;
klist = special[i];
knum = nspecial[i][0];
int molecular = atom->molecular;
if (molecular == 1) {
klist = atom->special[i];
knum = atom->nspecial[i][0];
} else {
if (atom->molindex[i] < 0) return 0.0;
int imol = atom->molindex[i];
int iatom = atom->molatom[i];
Molecule **onemols = atom->avec->onemols;
tlist = onemols[imol]->special[iatom];
knum = onemols[imol]->nspecial[iatom][0];
tagprev = atom->tag[i] - iatom - 1;
}
factor_hb = special_lj[sbmask(j)];
for (kk = 0; kk < knum; kk++) {
k = atom->map(klist[kk]);
if (molecular == 1) k = atom->map(klist[kk]);
else k = atom->map(tlist[kk]+tagprev);
if (k < 0) continue;
ktype = type[k];
m = type2param[itype][jtype][ktype];

48
src/MOLECULE/pair_hbond_dreiding_morse.cpp
View File

@ -21,6 +21,8 @@
#include "string.h"
#include "pair_hbond_dreiding_morse.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
@ -48,13 +50,15 @@ PairHbondDreidingMorse::PairHbondDreidingMorse(LAMMPS *lmp) :
void PairHbondDreidingMorse::compute(int eflag, int vflag)
{
int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype;
int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype,imol,iatom;
tagint tagprev;
double delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
double factor_hb,force_angle,force_kernel,force_switch,evdwl,ehbond;
double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
double fi[3],fj[3],delr1[3],delr2[3];
double r,dr,dexp,eng_morse,switch1,switch2;
int *ilist,*jlist,*numneigh,**firstneigh;
int *tlist;
tagint *klist;
evdwl = ehbond = 0.0;
@ -63,10 +67,15 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
tagint **special = atom->special;
int *type = atom->type;
int **nspecial = atom->nspecial;
int *type = atom->type;
double *special_lj = force->special_lj;
int molecular = atom->molecular;
Molecule **onemols = atom->avec->onemols;
inum = list->inum;
ilist = list->ilist;
@ -83,8 +92,17 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
i = ilist[ii];
itype = type[i];
if (!donor[itype]) continue;
if (molecular == 1) {
klist = special[i];
knum = nspecial[i][0];
} else {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
tlist = onemols[imol]->special[iatom];
knum = onemols[imol]->nspecial[iatom][0];
tagprev = tag[i] - iatom - 1;
}
jlist = firstneigh[i];
jnum = numneigh[i];
@ -102,7 +120,8 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
rsq = delx*delx + dely*dely + delz*delz;
for (kk = 0; kk < knum; kk++) {
k = atom->map(klist[kk]);
if (molecular == 1) k = atom->map(klist[kk]);
else k = atom->map(tlist[kk]+tagprev);
if (k < 0) continue;
ktype = type[k];
m = type2param[itype][jtype][ktype];
@ -357,16 +376,16 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype,
double &fforce)
{
int k,kk,ktype,knum,m;
tagint tagprev;
double eng,eng_morse,force_kernel,force_angle;
double rsq1,rsq2,r1,r2,c,s,ac,r,dr,dexp,factor_hb;
double switch1,switch2;
double delr1[3],delr2[3];
int *tlist;
tagint *klist;
double **x = atom->x;
tagint **special = atom->special;
int *type = atom->type;
int **nspecial = atom->nspecial;
double *special_lj = force->special_lj;
eng = 0.0;
@ -377,13 +396,26 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype,
if (!donor[itype]) return 0.0;
if (!acceptor[jtype]) return 0.0;
klist = special[i];
knum = nspecial[i][0];
int molecular = atom->molecular;
if (molecular == 1) {
klist = atom->special[i];
knum = atom->nspecial[i][0];
} else {
if (atom->molindex[i] < 0) return 0.0;
int imol = atom->molindex[i];
int iatom = atom->molatom[i];
Molecule **onemols = atom->avec->onemols;
tlist = onemols[imol]->special[iatom];
knum = onemols[imol]->nspecial[iatom][0];
tagprev = atom->tag[i] - iatom - 1;
}
factor_hb = special_lj[sbmask(j)];
for (kk = 0; kk < knum; kk++) {
k = atom->map(klist[kk]);
if (molecular == 1) k = atom->map(klist[kk]);
else k = atom->map(tlist[kk]+tagprev);
if (k < 0) continue;
ktype = type[k];
m = type2param[itype][jtype][ktype];

8
src/RIGID/fix_rigid_small.cpp
View File

@ -1288,7 +1288,7 @@ void FixRigidSmall::create_bodies()
// value = index into per-body data structure
// n = # of entries in hash
hash = new std::map<int,int>();
hash = new std::map<tagint,int>();
hash->clear();
// setup hash
@ -1296,7 +1296,7 @@ void FixRigidSmall::create_bodies()
// value = index into N-length data structure
// n = count of unique bodies my atoms are part of
int *molecule = atom->molecule;
tagint *molecule = atom->molecule;
n = 0;
for (i = 0; i < nlocal; i++) {
@ -1430,7 +1430,7 @@ void FixRigidSmall::create_bodies()
void FixRigidSmall::ring_bbox(int n, char *cbuf)
{
std::map<int,int> *hash = frsptr->hash;
std::map<tagint,int> *hash = frsptr->hash;
double **bbox = frsptr->bbox;
double *buf = (double *) cbuf;
@ -1462,7 +1462,7 @@ void FixRigidSmall::ring_bbox(int n, char *cbuf)
void FixRigidSmall::ring_nearest(int n, char *cbuf)
{
std::map<int,int> *hash = frsptr->hash;
std::map<tagint,int> *hash = frsptr->hash;
double **ctr = frsptr->ctr;
tagint *idclose = frsptr->idclose;
double *rsqclose = frsptr->rsqclose;

2
src/RIGID/fix_rigid_small.h
View File

@ -157,7 +157,7 @@ class FixRigidSmall : public Fix {
// class data used by ring communication callbacks
std::map<int,int> *hash;
std::map<tagint,int> *hash;
double **bbox;
double **ctr;
tagint *idclose;

222
src/RIGID/fix_shake.cpp
View File

@ -59,7 +59,8 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
// error check
if (atom->molecular == 0)
molecular = atom->molecular;
if (molecular == 0)
error->all(FLERR,"Cannot use fix shake with non-molecular system");
// perform initial allocation of atom-based arrays
@ -219,35 +220,23 @@ FixShake::~FixShake()
// set bond_type and angle_type back to positive for SHAKE clusters
// must set for all SHAKE bonds and angles stored by each atom
int **bond_type = atom->bond_type;
int **angle_type = atom->angle_type;
int nlocal = atom->nlocal;
int n;
for (int i = 0; i < nlocal; i++) {
if (shake_flag[i] == 0) continue;
else if (shake_flag[i] == 1) {
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = anglefind(i,shake_atom[i][1],shake_atom[i][2]);
if (n >= 0) angle_type[i][n] = -angle_type[i][n];
bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],1);
bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],1);
angletype_findset(i,shake_atom[i][1],shake_atom[i][2],1);
} else if (shake_flag[i] == 2) {
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],1);
} else if (shake_flag[i] == 3) {
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],1);
bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],1);
} else if (shake_flag[i] == 4) {
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = bondfind(i,shake_atom[i][0],shake_atom[i][3]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],1);
bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],1);
bondtype_findset(i,shake_atom[i][0],shake_atom[i][3],1);
}
}
@ -655,8 +644,9 @@ int FixShake::dof(int igroup)
void FixShake::find_clusters()
{
int i,j,m,n;
int i,j,m,n,imol,iatom;
int flag,flag_all,messtag,loop,nbuf,nbufmax,size;
tagint tagprev;
double massone;
tagint *buf;
MPI_Request request;
@ -664,19 +654,20 @@ void FixShake::find_clusters()
if (me == 0 && screen) fprintf(screen,"Finding SHAKE clusters ...\n");
// local copies of atom ptrs
onemols = atom->avec->onemols;
tagint *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
double *mass = atom->mass;
double *rmass = atom->rmass;
int **bond_type = atom->bond_type;
int **angle_type = atom->angle_type;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
int nlocal = atom->nlocal;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int nlocal = atom->nlocal;
int angles_allow = atom->avec->angles_allow;
// setup ring of procs
@ -701,7 +692,16 @@ void FixShake::find_clusters()
// -----------------------------------------------------
int max = 0;
if (molecular == 1) {
for (i = 0; i < nlocal; i++) max = MAX(max,nspecial[i][0]);
} else {
for (i = 0; i < nlocal; i++) {
imol = molindex[i];
if (imol < 0) continue;
iatom = molatom[i];
max = MAX(max,onemols[imol]->nspecial[iatom][0]);
}
}
int *npartner;
memory->create(npartner,nlocal,"shake:npartner");
@ -722,11 +722,23 @@ void FixShake::find_clusters()
// set npartner and partner_tag from special arrays
// -----------------------------------------------------
if (molecular == 1) {
for (i = 0; i < nlocal; i++) {
npartner[i] = nspecial[i][0];
for (j = 0; j < npartner[i]; j++)
partner_tag[i][j] = special[i][j];
}
} else {
for (i = 0; i < nlocal; i++) {
imol = molindex[i];
if (imol < 0) continue;
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
npartner[i] = onemols[imol]->nspecial[iatom][0];
for (j = 0; j < npartner[i]; j++)
partner_tag[i][j] = onemols[imol]->special[iatom][j] + tagprev;;
}
}
// -----------------------------------------------------
// set partner_mask, partner_type, partner_massflag, partner_bondtype
@ -758,11 +770,11 @@ void FixShake::find_clusters()
else massone = mass[type[m]];
partner_massflag[i][j] = masscheck(massone);
}
n = bondfind(i,tag[i],partner_tag[i][j]);
if (n >= 0) partner_bondtype[i][j] = bond_type[i][n];
n = bondtype_findset(i,tag[i],partner_tag[i][j],0);
if (n) partner_bondtype[i][j] = n;
else {
n = bondfind(m,tag[i],partner_tag[i][j]);
if (n >= 0) partner_bondtype[i][j] = bond_type[m][n];
n = bondtype_findset(m,tag[i],partner_tag[i][j],0);
if (n) partner_bondtype[i][j] = n;
}
} else nbuf += nper;
}
@ -782,8 +794,8 @@ void FixShake::find_clusters()
buf[size+2] = 0;
buf[size+3] = 0;
buf[size+4] = 0;
n = bondfind(i,tag[i],partner_tag[i][j]);
if (n >= 0) buf[size+5] = bond_type[i][n];
n = bondtype_findset(i,tag[i],partner_tag[i][j],0);
if (n) buf[size+5] = n;
else buf[size+5] = 0;
size += nper;
}
@ -1006,12 +1018,11 @@ void FixShake::find_clusters()
}
if (nshake[i] == 2 && angles_allow) {
n = anglefind(i,shake_atom[i][1],shake_atom[i][2]);
if (n < 0) continue;
if (angle_type[i][n] < 0) continue;
if (angle_flag[angle_type[i][n]]) {
n = angletype_findset(i,shake_atom[i][1],shake_atom[i][2],0);
if (n <= 0) continue;
if (angle_flag[n]) {
shake_flag[i] = 1;
shake_type[i][2] = angle_type[i][n];
shake_type[i][2] = n;
}
}
}
@ -1099,27 +1110,18 @@ void FixShake::find_clusters()
for (i = 0; i < nlocal; i++) {
if (shake_flag[i] == 0) continue;
else if (shake_flag[i] == 1) {
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = anglefind(i,shake_atom[i][1],shake_atom[i][2]);
if (n >= 0) angle_type[i][n] = -angle_type[i][n];
bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1);
angletype_findset(i,shake_atom[i][1],shake_atom[i][2],-1);
} else if (shake_flag[i] == 2) {
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
} else if (shake_flag[i] == 3) {
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1);
} else if (shake_flag[i] == 4) {
n = bondfind(i,shake_atom[i][0],shake_atom[i][1]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = bondfind(i,shake_atom[i][0],shake_atom[i][2]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
n = bondfind(i,shake_atom[i][0],shake_atom[i][3]);
if (n >= 0) bond_type[i][n] = -bond_type[i][n];
bondtype_findset(i,shake_atom[i][0],shake_atom[i][1],-1);
bondtype_findset(i,shake_atom[i][0],shake_atom[i][2],-1);
bondtype_findset(i,shake_atom[i][0],shake_atom[i][3],-1);
}
}
@ -1175,7 +1177,6 @@ void FixShake::ring_bonds(int ndatum, char *cbuf)
double *rmass = atom->rmass;
double *mass = atom->mass;
int *mask = atom->mask;
int **bond_type = atom->bond_type;
int *type = atom->type;
int nlocal = atom->nlocal;
int nmass = fsptr->nmass;
@ -1195,8 +1196,8 @@ void FixShake::ring_bonds(int ndatum, char *cbuf)
buf[i+4] = fsptr->masscheck(massone);
}
if (buf[i+5] == 0) {
n = fsptr->bondfind(m,buf[i],buf[i+1]);
if (n >= 0) buf[i+5] = bond_type[m][n];
n = fsptr->bondtype_findset(m,buf[i],buf[i+1],0);
if (n) buf[i+5] = n;
}
}
}
@ -2236,45 +2237,116 @@ void FixShake::stats()
}
/* ----------------------------------------------------------------------
find a bond between global tags n1 and n2 stored with local atom i
return -1 if don't find it
return bond index if do find it
find a bond between global atom IDs n1 and n2 stored with local atom i
if find it:
if setflag = 0, return bond type
if setflag = -1/1, set bond type to negative/positive and return 0
if do not find it, return 0
------------------------------------------------------------------------- */
int FixShake::bondfind(int i, tagint n1, tagint n2)
int FixShake::bondtype_findset(int i, tagint n1, tagint n2, int setflag)
{
int m,nbonds;
int *btype;
if (molecular == 1) {
tagint *tag = atom->tag;
tagint **bond_atom = atom->bond_atom;
int nbonds = atom->num_bond[i];
nbonds = atom->num_bond[i];
int m;
for (m = 0; m < nbonds; m++) {
if (n1 == tag[i] && n2 == bond_atom[i][m]) break;
if (n1 == bond_atom[i][m] && n2 == tag[i]) break;
}
if (m < nbonds) return m;
return -1;
} else {
int imol = atom->molindex[i];
int iatom = atom->molatom[i];
tagint *tag = atom->tag;
tagint tagprev = tag[i] - iatom - 1;
int *batom = onemols[imol]->bond_atom[iatom];
btype = onemols[imol]->bond_type[iatom];
nbonds = onemols[imol]->num_bond[iatom];
for (m = 0; m < nbonds; m++) {
if (n1 == tag[i] && n2 == batom[m]+tagprev) break;
if (n1 == batom[m]+tagprev && n2 == tag[i]) break;
}
}
if (m < nbonds) {
if (setflag == 0) {
if (molecular == 1) return atom->bond_type[i][m];
else return btype[m];
}
if (molecular == 1) {
if ((setflag < 0 && atom->bond_type[i][m] > 0) ||
(setflag > 0 && atom->bond_type[i][m] < 0))
atom->bond_type[i][m] = -atom->bond_type[i][m];
} else {
if ((setflag < 0 && btype[m] > 0) ||
(setflag > 0 && btype[m] < 0)) btype[m] = -btype[m];
}
}
return 0;
}
/* ----------------------------------------------------------------------
find an angle with global end atoms n1 and n2 stored with local atom i
return -1 if don't find it
return angle index if do find it
find an angle with global end atom IDs n1 and n2 stored with local atom i
if find it:
if setflag = 0, return angle type
if setflag = -1/1, set angle type to negative/positive and return 0
if do not find it, return 0
------------------------------------------------------------------------- */
int FixShake::anglefind(int i, tagint n1, tagint n2)
int FixShake::angletype_findset(int i, tagint n1, tagint n2, int setflag)
{
int m,nangles;
int *atype;
if (molecular == 1) {
tagint **angle_atom1 = atom->angle_atom1;
tagint **angle_atom3 = atom->angle_atom3;
int nangles = atom->num_angle[i];
nangles = atom->num_angle[i];
int m;
for (m = 0; m < nangles; m++) {
if (n1 == angle_atom1[i][m] && n2 == angle_atom3[i][m]) break;
if (n1 == angle_atom3[i][m] && n2 == angle_atom1[i][m]) break;
}
if (m < nangles) return m;
return -1;
} else {
int imol = atom->molindex[i];
int iatom = atom->molatom[i];
tagint *tag = atom->tag;
tagint tagprev = tag[i] - iatom - 1;
int *aatom1 = onemols[imol]->angle_atom1[iatom];
int *aatom3 = onemols[imol]->angle_atom3[iatom];
atype = onemols[imol]->angle_type[iatom];
nangles = onemols[imol]->num_angle[iatom];
for (m = 0; m < nangles; m++) {
if (n1 == aatom1[m]+tagprev && n2 == aatom3[m]+tagprev) break;
if (n1 == aatom3[m]+tagprev && n2 == aatom1[m]+tagprev) break;
}
}
if (m < nangles) {
if (setflag == 0) {
if (molecular == 1) return atom->angle_type[i][m];
else return atype[m];
}
if (molecular == 1) {
if ((setflag < 0 && atom->angle_type[i][m] > 0) ||
(setflag > 0 && atom->angle_type[i][m] < 0))
atom->angle_type[i][m] = -atom->angle_type[i][m];
} else {
if ((setflag < 0 && atype[m] > 0) ||
(setflag > 0 && atype[m] < 0)) atype[m] = -atype[m];
}
}
return 0;
}
/* ----------------------------------------------------------------------

10
src/RIGID/fix_shake.h
View File

@ -64,6 +64,7 @@ class FixShake : public Fix {
double *mass_list; // constrain bonds to these masses
int nmass; // # of masses in mass_list
int molecular; // copy of atom->molecular
double *bond_distance,*angle_distance; // constraint distances
int ifix_respa; // rRESPA fix needed by SHAKE
@ -102,9 +103,8 @@ class FixShake : public Fix {
double *a_ave,*a_max,*a_min;
double *a_ave_all,*a_max_all,*a_min_all;
// molecules added on-the-fly with SHAKE constraints
class Molecule *onemol;
class Molecule **onemols; // atom style template pointer
class Molecule *onemol; // molecule added on-the-fly
void find_clusters();
int masscheck(double);
@ -115,8 +115,8 @@ class FixShake : public Fix {
void shake4(int);
void shake3angle(int);
void stats();
int bondfind(int, tagint, tagint);
int anglefind(int, tagint, tagint);
int bondtype_findset(int, tagint, tagint, int);
int angletype_findset(int, tagint, tagint, int);
// static variable for ring communication callback to access class data
// callback functions for ring communication

42
src/atom.cpp
View File

@ -74,6 +74,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
x = v = f = NULL;
molecule = NULL;
molindex = molatom = NULL;
q = NULL;
mu = NULL;
omega = angmom = torque = NULL;
@ -215,6 +216,8 @@ Atom::~Atom()
memory->destroy(erforce);
memory->destroy(molecule);
memory->destroy(molindex);
memory->destroy(molatom);
memory->destroy(nspecial);
memory->destroy(special);
@ -321,7 +324,8 @@ void Atom::create_avec(const char *style, int narg, char **arg, char *suffix)
int sflag;
avec = new_avec(style,suffix,sflag);
avec->settings(narg,arg);
avec->store_args(narg,arg);
avec->process_args(narg,arg);
avec->grow(1);
nmax = 0;
avec->reset();
@ -1336,30 +1340,33 @@ int Atom::find_molecule(char *id)
void Atom::add_molecule_atom(Molecule *onemol, int iatom,
int ilocal, tagint offset)
{
if (onemol->qflag) q[ilocal] = onemol->q[iatom];
if (onemol->radiusflag) radius[ilocal] = onemol->radius[iatom];
if (onemol->rmassflag) rmass[ilocal] = onemol->rmass[iatom];
if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom];
if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom];
if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom];
else if (rmass_flag)
rmass[ilocal] = 4.0*MY_PI/3.0 *
radius[ilocal]*radius[ilocal]*radius[ilocal];
if (onemol->bondflag) {
num_bond[ilocal] = onemol->num_bond[iatom];
if (molecular != 1) return;
// add bond topology info
// for molecular atom styles, but not atom style template
if (avec->bonds_allow) num_bond[ilocal] = onemol->num_bond[iatom];
for (int i = 0; i < num_bond[ilocal]; i++) {
bond_type[ilocal][i] = onemol->bond_type[iatom][i];
bond_atom[ilocal][i] = onemol->bond_atom[iatom][i] + offset;
}
}
if (onemol->angleflag) {
num_angle[ilocal] = onemol->num_angle[iatom];
if (avec->angles_allow) num_angle[ilocal] = onemol->num_angle[iatom];
for (int i = 0; i < num_angle[ilocal]; i++) {
angle_type[ilocal][i] = onemol->angle_type[iatom][i];
angle_atom1[ilocal][i] = onemol->angle_atom1[iatom][i] + offset;
angle_atom2[ilocal][i] = onemol->angle_atom2[iatom][i] + offset;
angle_atom3[ilocal][i] = onemol->angle_atom3[iatom][i] + offset;
}
}
if (onemol->dihedralflag) {
if (avec->dihedrals_allow)
num_dihedral[ilocal] = onemol->num_dihedral[iatom];
for (int i = 0; i < num_dihedral[ilocal]; i++) {
dihedral_type[ilocal][i] = onemol->dihedral_type[iatom][i];
@ -1368,8 +1375,8 @@ void Atom::add_molecule_atom(Molecule *onemol, int iatom,
dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset;
dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset;
}
}
if (onemol->improperflag) {
if (avec->impropers_allow)
num_improper[ilocal] = onemol->num_improper[iatom];
for (int i = 0; i < num_improper[ilocal]; i++) {
improper_type[ilocal][i] = onemol->improper_type[iatom][i];
@ -1378,17 +1385,8 @@ void Atom::add_molecule_atom(Molecule *onemol, int iatom,
improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset;
improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset;
}
}
// error check against maxspecial in case user has not done one of these:
// create_box extra/special/per/atom N
// read_data extra special per atom N
// special_bonds extra N
// if explicitly used special_bonds, may not have maintained extra
if (onemol->specialflag) {
if (onemol->maxspecial > maxspecial)
error->one(FLERR,"Molecule file special bond counts are too large");
nspecial[ilocal][0] = onemol->nspecial[iatom][0];
nspecial[ilocal][1] = onemol->nspecial[iatom][1];
int n = nspecial[ilocal][2] = onemol->nspecial[iatom][2];

3
src/atom.h
View File

@ -30,7 +30,8 @@ class Atom : protected Pointers {
int nlocal,nghost; // # of owned and ghost atoms on this proc
int nmax; // max # of owned+ghost in arrays on this proc
int tag_enable; // 0/1 if atom ID tags are defined
int molecular; // 0 = atomic, 1 = molecular system
int molecular; // 0 = atomic, 1 = standard molecular system,
// 2 = molecule template system
bigint nbonds,nangles,ndihedrals,nimpropers;
int ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;

29
src/atom_vec.cpp
View File

@ -11,6 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "atom_vec.h"
#include "atom.h"
@ -29,13 +30,39 @@ AtomVec::AtomVec(LAMMPS *lmp) : Pointers(lmp)
mass_type = dipole_type = 0;
size_data_bonus = 0;
cudable = false;
nargcopy = 0;
argcopy = NULL;
}
/* ---------------------------------------------------------------------- */
AtomVec::~AtomVec()
{
for (int i = 0; i < nargcopy; i++) delete [] argcopy[i];
delete [] argcopy;
}
/* ----------------------------------------------------------------------
make copy of args for use by restart & replicate
------------------------------------------------------------------------- */
void AtomVec::store_args(int narg, char **arg)
{
nargcopy = narg;
argcopy = new char*[nargcopy];
for (int i = 0; i < nargcopy; i++) {
int n = strlen(arg[i]) + 1;
argcopy[i] = new char[n];
strcpy(argcopy[i],arg[i]);
}
}
/* ----------------------------------------------------------------------
no additional args by default
------------------------------------------------------------------------- */
void AtomVec::settings(int narg, char **arg)
void AtomVec::process_args(int narg, char **arg)
{
if (narg) error->all(FLERR,"Invalid atom_style command");
}

14
src/atom_vec.h
View File

@ -39,12 +39,19 @@ class AtomVec : protected Pointers {
int size_data_bonus; // number of values in Bonus line
int xcol_data; // column (1-N) where x is in Atom line
class Molecule **onemols; // list of molecules for style template
int nset; // # of molecules in list
int cudable; // 1 if atom style is CUDA-enabled
int *maxsend; // CUDA-specific variable
int nargcopy; // copy of command-line args for atom_style command
char **argcopy; // used when AtomVec is realloced (restart,replicate)
AtomVec(class LAMMPS *);
virtual ~AtomVec() {}
virtual void settings(int, char **);
virtual ~AtomVec();
void store_args(int, char **);
virtual void process_args(int, char **);
virtual void init();
virtual void grow(int) = 0;
@ -78,9 +85,6 @@ class AtomVec : protected Pointers {
virtual int pack_restart(int, double *) = 0;
virtual int unpack_restart(double *) = 0;
virtual void write_restart_settings(FILE *) {}
virtual void read_restart_settings(FILE *) {}
virtual void create_atom(int, double *) = 0;
virtual void data_atom(double *, imageint, char **) = 0;

57
src/atom_vec_body.cpp
View File

@ -57,10 +57,6 @@ AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp)
bptr = NULL;
nargcopy = 0;
argcopy = NULL;
copyflag = 1;
if (sizeof(double) == sizeof(int)) intdoubleratio = 1;
else if (sizeof(double) == 2*sizeof(int)) intdoubleratio = 2;
else error->all(FLERR,"Internal error in atom_style body");
@ -78,9 +74,6 @@ AtomVecBody::~AtomVecBody()
memory->sfree(bonus);
delete bptr;
for (int i = 0; i < nargcopy; i++) delete [] argcopy[i];
delete [] argcopy;
}
/* ----------------------------------------------------------------------
@ -89,7 +82,7 @@ AtomVecBody::~AtomVecBody()
set size_forward and size_border to max sizes
------------------------------------------------------------------------- */
void AtomVecBody::settings(int narg, char **arg)
void AtomVecBody::process_args(int narg, char **arg)
{
if (narg < 1) error->all(FLERR,"Invalid atom_style body command");
@ -114,19 +107,6 @@ void AtomVecBody::settings(int narg, char **arg)
size_forward = 7 + bptr->size_forward;
size_border = 16 + bptr->size_border;
// make copy of args if called externally, so can write to restart file
// make no copy of args if called from read_restart()
if (copyflag) {
nargcopy = narg;
argcopy = new char*[nargcopy];
for (int i = 0; i < nargcopy; i++) {
int n = strlen(arg[i]) + 1;
argcopy[i] = new char[n];
strcpy(argcopy[i],arg[i]);
}
}
}
/* ----------------------------------------------------------------------
@ -1229,41 +1209,6 @@ int AtomVecBody::unpack_restart(double *buf)
return m;
}
/* ---------------------------------------------------------------------- */
void AtomVecBody::write_restart_settings(FILE *fp)
{
fwrite(&nargcopy,sizeof(int),1,fp);
for (int i = 0; i < nargcopy; i++) {
int n = strlen(argcopy[i]) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(argcopy[i],sizeof(char),n,fp);
}
}
/* ---------------------------------------------------------------------- */
void AtomVecBody::read_restart_settings(FILE *fp)
{
int n;
int me = comm->me;
if (me == 0) fread(&nargcopy,sizeof(int),1,fp);
MPI_Bcast(&nargcopy,1,MPI_INT,0,world);
argcopy = new char*[nargcopy];
for (int i = 0; i < nargcopy; i++) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
argcopy[i] = new char[n];
if (me == 0) fread(argcopy[i],sizeof(char),n,fp);
MPI_Bcast(argcopy[i],n,MPI_CHAR,0,world);
}
copyflag = 0;
settings(nargcopy,argcopy);
}
/* ----------------------------------------------------------------------
create one atom of itype at coord
set other values to defaults

8
src/atom_vec_body.h
View File

@ -42,7 +42,7 @@ class AtomVecBody : public AtomVec {
AtomVecBody(class LAMMPS *);
~AtomVecBody();
void settings(int, char **);
void process_args(int, char **);
void grow(int);
void grow_reset();
void copy(int, int, int);
@ -67,8 +67,6 @@ class AtomVecBody : public AtomVec {
int size_restart();
int pack_restart(int, double *);
int unpack_restart(double *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void create_atom(int, double *);
void data_atom(double *, imageint, char **);
int data_atom_hybrid(int, char **);
@ -99,10 +97,6 @@ class AtomVecBody : public AtomVec {
int *body;
int nlocal_bonus,nghost_bonus,nmax_bonus;
int nargcopy; // copy of command-line args
char **argcopy; // for writing to restart file
int copyflag;
int intdoubleratio; // sizeof(double) / sizeof(int)
MyPoolChunk<int> *icp;

27
src/atom_vec_hybrid.cpp
View File

@ -43,7 +43,7 @@ AtomVecHybrid::~AtomVecHybrid()
process sub-style args
------------------------------------------------------------------------- */
void AtomVecHybrid::settings(int narg, char **arg)
void AtomVecHybrid::process_args(int narg, char **arg)
{
// build list of all known atom styles
@ -55,7 +55,7 @@ void AtomVecHybrid::settings(int narg, char **arg)
keywords = new char*[narg];
// allocate each sub-style
// call settings() with set of args that are not atom style names
// call process_args() with set of args that are not atom style names
// use known_style() to determine which args these are
int i,jarg,dummy;
@ -73,7 +73,7 @@ void AtomVecHybrid::settings(int narg, char **arg)
strcpy(keywords[nstyles],arg[iarg]);
jarg = iarg + 1;
while (jarg < narg && !known_style(arg[jarg])) jarg++;
styles[nstyles]->settings(jarg-iarg-1,&arg[iarg+1]);
styles[nstyles]->process_args(jarg-iarg-1,&arg[iarg+1]);
iarg = jarg;
nstyles++;
}
@ -97,7 +97,12 @@ void AtomVecHybrid::settings(int narg, char **arg)
xcol_data = 3;
for (int k = 0; k < nstyles; k++) {
if ((styles[k]->molecular == 1 && molecular == 2) ||
(styles[k]->molecular == 2 && molecular == 1))
error->all(FLERR,"Cannot mix molecular and molecule template "
"atom styles");
molecular = MAX(molecular,styles[k]->molecular);
bonds_allow = MAX(bonds_allow,styles[k]->bonds_allow);
angles_allow = MAX(angles_allow,styles[k]->angles_allow);
dihedrals_allow = MAX(dihedrals_allow,styles[k]->dihedrals_allow);
@ -820,22 +825,6 @@ int AtomVecHybrid::unpack_restart(double *buf)
return m;
}
/* ---------------------------------------------------------------------- */
void AtomVecHybrid::write_restart_settings(FILE *fp)
{
for (int k = 0; k < nstyles; k++)
styles[k]->write_restart_settings(fp);
}
/* ---------------------------------------------------------------------- */
void AtomVecHybrid::read_restart_settings(FILE *fp)
{
for (int k = 0; k < nstyles; k++)
styles[k]->read_restart_settings(fp);
}
/* ----------------------------------------------------------------------
create one atom of itype at coord
create each sub-style one after the other

4
src/atom_vec_hybrid.h
View File

@ -33,7 +33,7 @@ class AtomVecHybrid : public AtomVec {
AtomVecHybrid(class LAMMPS *);
~AtomVecHybrid();
void settings(int, char **);
void process_args(int, char **);
void init();
void grow(int);
void grow_reset();
@ -54,8 +54,6 @@ class AtomVecHybrid : public AtomVec {
int size_restart();
int pack_restart(int, double *);
int unpack_restart(double *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void create_atom(int, double *);
void data_atom(double *, imageint, char **);
int data_atom_hybrid(int, char **) {return 0;}

2
src/atom_vec_line.h
View File

@ -88,7 +88,7 @@ class AtomVecLine : public AtomVec {
int *type,*mask;
imageint *image;
double **x,**v,**f;
int *molecule;
tagint *molecule;
double *rmass;
double **omega,**torque;
int *line;

2
src/atom_vec_tri.h
View File

@ -90,7 +90,7 @@ class AtomVecTri : public AtomVec {
int *type,*mask;
imageint *image;
double **x,**v,**f;
int *molecule;
tagint *molecule;
double *rmass;
double **angmom,**torque;
int *tri;

43
src/compute_angle_local.cpp
View File

@ -16,6 +16,7 @@
#include "compute_angle_local.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h"
#include "force.h"
@ -98,7 +99,7 @@ void ComputeAngleLocal::compute_local()
/* ----------------------------------------------------------------------
count angles and compute angle info on this proc
only count angle once if newton_angle is off
only count if 2nd atom is the one storing the angle
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
@ -109,20 +110,27 @@ void ComputeAngleLocal::compute_local()
int ComputeAngleLocal::compute_angles(int flag)
{
int i,m,n,atom1,atom2,atom3;
int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype;
tagint tagprev;
double delx1,dely1,delz1,delx2,dely2,delz2;
double rsq1,rsq2,r1,r2,c;
double *tbuf,*ebuf;
double **x = atom->x;
tagint *tag = atom->tag;
int *num_angle = atom->num_angle;
tagint **angle_atom1 = atom->angle_atom1;
tagint **angle_atom2 = atom->angle_atom2;
tagint **angle_atom3 = atom->angle_atom3;
int **angle_type = atom->angle_type;
tagint *tag = atom->tag;
int *mask = atom->mask;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (flag) {
if (nvalues == 1) {
@ -141,13 +149,32 @@ int ComputeAngleLocal::compute_angles(int flag)
m = n = 0;
for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue;
for (i = 0; i < num_angle[atom2]; i++) {
if (molecular == 1) na = num_angle[atom2];
else {
if (molindex[atom2] < 0) continue;
imol = molindex[atom2];
iatom = molatom[atom2];
na = onemols[imol]->num_angle[iatom];
}
for (i = 0; i < na; i++) {
if (molecular == 1) {
if (tag[atom2] != angle_atom2[atom2][i]) continue;
atype = angle_type[atom2][i];
atom1 = atom->map(angle_atom1[atom2][i]);
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
atom3 = atom->map(angle_atom3[atom2][i]);
} else {
if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue;
tagprev = tag[atom1] - iatom - 1;
atype = atom->map(onemols[imol]->angle_type[atom2][i]);
atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i]+tagprev);
atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i]+tagprev);
}
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
if (angle_type[atom2][i] == 0) continue;
if (atype == 0) continue;
if (flag) {
if (tflag >= 0) {
@ -177,8 +204,8 @@ int ComputeAngleLocal::compute_angles(int flag)
}
if (eflag >= 0) {
if (angle_type[atom2][i] > 0)
ebuf[n] = angle->single(angle_type[atom2][i],atom1,atom2,atom3);
if (atype > 0)
ebuf[n] = angle->single(atype,atom1,atom2,atom3);
else ebuf[n] = 0.0;
}
n += nvalues;

37
src/compute_bond_local.cpp
View File

@ -16,6 +16,7 @@
#include "compute_bond_local.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h"
#include "force.h"
@ -115,19 +116,26 @@ void ComputeBondLocal::compute_local()
int ComputeBondLocal::compute_bonds(int flag)
{
int i,m,n,atom1,atom2;
int i,m,n,nb,atom1,atom2,imol,iatom,btype;
tagint tagprev;
double delx,dely,delz,rsq;
double *dbuf,*ebuf,*fbuf;
double *ptr;
double **x = atom->x;
tagint *tag = atom->tag;
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
tagint *tag = atom->tag;
int *mask = atom->mask;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
int molecular = atom->molecular;
Bond *bond = force->bond;
double eng,fbond;
@ -135,11 +143,28 @@ int ComputeBondLocal::compute_bonds(int flag)
m = n = 0;
for (atom1 = 0; atom1 < nlocal; atom1++) {
if (!(mask[atom1] & groupbit)) continue;
for (i = 0; i < num_bond[atom1]; i++) {
if (molecular == 1) nb = num_bond[atom1];
else {
if (molindex[atom1] < 0) continue;
imol = molindex[atom1];
iatom = molatom[atom1];
nb = onemols[imol]->num_bond[iatom];
}
for (i = 0; i < nb; i++) {
if (molecular == 1) {
btype = bond_type[atom1][i];
atom2 = atom->map(bond_atom[atom1][i]);
} else {
tagprev = tag[atom1] - iatom - 1;
btype = atom->map(onemols[imol]->bond_type[atom1][i]);
atom2 = atom->map(onemols[imol]->bond_atom[atom1][i]+tagprev);
}
if (atom2 < 0 || !(mask[atom2] & groupbit)) continue;
if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
if (bond_type[atom1][i] == 0) continue;
if (btype == 0) continue;
if (flag) {
delx = x[atom1][0] - x[atom2][0];
@ -149,8 +174,8 @@ int ComputeBondLocal::compute_bonds(int flag)
rsq = delx*delx + dely*dely + delz*delz;
if (singleflag) {
if (bond_type[atom1][i] > 0)
eng = bond->single(bond_type[atom1][i],rsq,atom1,atom2,fbond);
if (btype > 0)
eng = bond->single(btype,rsq,atom1,atom2,fbond);
else eng = fbond = 0.0;
}

36
src/compute_dihedral_local.cpp
View File

@ -16,6 +16,7 @@
#include "compute_dihedral_local.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h"
#include "force.h"
@ -107,21 +108,28 @@ void ComputeDihedralLocal::compute_local()
int ComputeDihedralLocal::compute_dihedrals(int flag)
{
int i,m,n,atom1,atom2,atom3,atom4;
int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom;
tagint tagprev;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
double s,c;
double *pbuf;
double **x = atom->x;
tagint *tag = atom->tag;
int *num_dihedral = atom->num_dihedral;
tagint **dihedral_atom1 = atom->dihedral_atom1;
tagint **dihedral_atom2 = atom->dihedral_atom2;
tagint **dihedral_atom3 = atom->dihedral_atom3;
tagint **dihedral_atom4 = atom->dihedral_atom4;
tagint *tag = atom->tag;
int *mask = atom->mask;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (flag) {
if (nvalues == 1) {
@ -135,13 +143,31 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
m = n = 0;
for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue;
for (i = 0; i < num_dihedral[atom2]; i++) {
if (molecular == 1) nd = num_dihedral[atom2];
else {
if (molindex[atom2] < 0) continue;
imol = molindex[atom2];
iatom = molatom[atom2];
nd = onemols[imol]->num_dihedral[iatom];
}
for (i = 0; i < nd; i++) {
if (molecular == 1) {
if (tag[atom2] != dihedral_atom2[atom2][i]) continue;
atom1 = atom->map(dihedral_atom1[atom2][i]);
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
atom3 = atom->map(dihedral_atom3[atom2][i]);
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
atom4 = atom->map(dihedral_atom4[atom2][i]);
} else {
if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue;
tagprev = tag[atom1] - iatom - 1;
atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i]+tagprev);
atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i]+tagprev);
atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i]+tagprev);
}
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
if (atom4 < 0 || !(mask[atom4] & groupbit)) continue;
if (flag) {

36
src/compute_improper_local.cpp
View File

@ -16,6 +16,7 @@
#include "compute_improper_local.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h"
#include "force.h"
@ -108,21 +109,28 @@ void ComputeImproperLocal::compute_local()
int ComputeImproperLocal::compute_impropers(int flag)
{
int i,m,n,atom1,atom2,atom3,atom4;
int i,m,n,ni,atom1,atom2,atom3,atom4,imol,iatom;
tagint tagprev;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z;
double ss1,ss2,ss3,r1,r2,r3,c0,c1,c2,s1,s2;
double s12,c;
double *cbuf;
double **x = atom->x;
tagint *tag = atom->tag;
int *num_improper = atom->num_improper;
tagint **improper_atom1 = atom->improper_atom1;
tagint **improper_atom2 = atom->improper_atom2;
tagint **improper_atom3 = atom->improper_atom3;
tagint **improper_atom4 = atom->improper_atom4;
tagint *tag = atom->tag;
int *mask = atom->mask;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (flag) {
if (nvalues == 1) {
@ -136,13 +144,31 @@ int ComputeImproperLocal::compute_impropers(int flag)
m = n = 0;
for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue;
for (i = 0; i < num_improper[atom2]; i++) {
if (molecular == 1) ni = num_improper[atom2];
else {
if (molindex[atom2] < 0) continue;
imol = molindex[atom2];
iatom = molatom[atom2];
ni = onemols[imol]->num_improper[iatom];
}
for (i = 0; i < ni; i++) {
if (molecular == 1) {
if (tag[atom2] != improper_atom2[atom2][i]) continue;
atom1 = atom->map(improper_atom1[atom2][i]);
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
atom3 = atom->map(improper_atom3[atom2][i]);
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
atom4 = atom->map(improper_atom4[atom2][i]);
} else {
if (tag[atom2] != onemols[imol]->improper_atom2[atom2][i]) continue;
tagprev = tag[atom1] - iatom - 1;
atom1 = atom->map(onemols[imol]->improper_atom1[atom2][i]+tagprev);
atom3 = atom->map(onemols[imol]->improper_atom3[atom2][i]+tagprev);
atom4 = atom->map(onemols[imol]->improper_atom4[atom2][i]+tagprev);
}
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
if (atom4 < 0 || !(mask[atom4] & groupbit)) continue;
if (flag) {

2
src/compute_property_atom.cpp
View File

@ -449,7 +449,7 @@ void ComputePropertyAtom::pack_id(int n)
void ComputePropertyAtom::pack_molecule(int n)
{
int *molecule = atom->molecule;
tagint *molecule = atom->molecule;
int *mask = atom->mask;
int nlocal = atom->nlocal;

5
src/compute_property_local.cpp
View File

@ -212,6 +212,11 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) :
// error check
if (atom->molecular == 2 && (kindflag == BOND || kindflag == ANGLE ||
kindflag == DIHEDRAL || kindflag == IMPROPER))
error->all(FLERR,"Compute property/local does not (yet) work "
"with atom_style template");
if (kindflag == BOND && atom->avec->bonds_allow == 0)
error->all(FLERR,
"Compute property/local for property that isn't allocated");

29
src/create_atoms.cpp
View File

@ -160,11 +160,12 @@ void CreateAtoms::command(int narg, char **arg)
error->all(FLERR,"Create_atoms molecule must have coordinates");
if (onemol->typeflag == 0)
error->all(FLERR,"Create_atoms molecule must have atom types");
if (ntype+onemol->maxtype <= 0 || ntype+onemol->maxtype > atom->ntypes)
if (ntype+onemol->ntypes <= 0 || ntype+onemol->ntypes > atom->ntypes)
error->all(FLERR,"Invalid atom type in create_atoms mol command");
if (onemol->tag_require && !atom->tag_enable)
error->all(FLERR,
"Create_atoms molecule has atom IDs, but system does not");
onemol->check_attributes(0);
// create_atoms uses geoemetric center of molecule for insertion
@ -291,6 +292,19 @@ void CreateAtoms::command(int narg, char **arg)
int molcreate = (atom->nlocal - nlocal_previous) / onemol->natoms;
// for atom style template systems, increment total bonds,angles,etc
if (atom->molecular == 2) {
bigint nmolme = molcreate;
bigint nmoltotal;
MPI_Allreduce(&nmolme,&nmoltotal,1,MPI_LMP_BIGINT,MPI_SUM,world);
atom->nbonds += nmoltotal * onemol->nbonds;
atom->nangles += nmoltotal * onemol->nangles;
atom->ndihedrals += nmoltotal * onemol->ndihedrals;
atom->nimpropers += nmoltotal * onemol->nimpropers;
}
// if atom style template
// maxmol = max molecule ID across all procs, for previous atoms
// moloffset = max molecule ID for all molecules owned by previous procs
// including molecules existing before this creation
@ -332,12 +346,17 @@ void CreateAtoms::command(int narg, char **arg)
tagint **improper_atom4 = atom->improper_atom4;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
int molecular = atom->molecular;
int ilocal = nlocal_previous;
for (int i = 0; i < molcreate; i++) {
if (tag) offset = tag[ilocal]-1;
for (int m = 0; m < natoms; m++) {
if (molecule) molecule[ilocal] = moloffset + i+1;
if (molecular) molecule[ilocal] = moloffset + i+1;
if (molecular == 2) {
atom->molindex[ilocal] = 0;
atom->molatom[ilocal] = m;
} else if (molecular) {
if (onemol->bondflag)
for (int j = 0; j < num_bond[ilocal]; j++)
bond_atom[ilocal][j] += offset;
@ -364,6 +383,7 @@ void CreateAtoms::command(int narg, char **arg)
if (onemol->specialflag)
for (int j = 0; j < nspecial[ilocal][2]; j++)
special[ilocal][j] += offset;
}
ilocal++;
}
}
@ -400,11 +420,12 @@ void CreateAtoms::command(int narg, char **arg)
}
// for MOLECULE mode:
// create special bond lists for molecular systems
// create special bond lists for molecular systems,
// but not for atom style template
// only if onemol added bonds but not special info
if (mode == MOLECULE) {
if (atom->molecular && onemol->bondflag && !onemol->specialflag) {
if (atom->molecular == 1 && onemol->bondflag && !onemol->specialflag) {
Special special(lmp);
special.build();
}

4
src/delete_atoms.cpp
View File

@ -184,7 +184,7 @@ void DeleteAtoms::delete_region(int narg, char **arg)
memory->create(list,n,"delete_atoms:list");
n = 0;
std::map<int,int>::iterator pos;
std::map<tagint,int>::iterator pos;
for (pos = hash->begin(); pos != hash->end(); ++pos) list[n++] = pos->first;
cptr = this;
@ -202,7 +202,7 @@ void DeleteAtoms::delete_region(int narg, char **arg)
void DeleteAtoms::molring(int n, char *cbuf)
{
tagint *list = (int *) cbuf;
tagint *list = (tagint *) cbuf;
int *dlist = cptr->dlist;
std::map<tagint,int> *hash = cptr->hash;
int nlocal = cptr->atom->nlocal;

2
src/delete_atoms.h
View File

@ -33,7 +33,7 @@ class DeleteAtoms : protected Pointers {
private:
int *dlist;
int compress_flag,mol_flag;
std::map<int,int> *hash;
std::map<tagint,int> *hash;
void delete_group(int, char **);
void delete_region(int, char **);

3
src/delete_bonds.cpp
View File

@ -41,8 +41,9 @@ void DeleteBonds::command(int narg, char **arg)
error->all(FLERR,"Delete_bonds command before simulation box is defined");
if (atom->natoms == 0)
error->all(FLERR,"Delete_bonds command with no atoms existing");
if (atom->molecular == 0)
if (atom->molecular != 1)
error->all(FLERR,"Cannot use delete_bonds with non-molecular system");
if (narg < 2) error->all(FLERR,"Illegal delete_bonds command");
// init entire system since comm->borders is done

65
src/domain.cpp
View File

@ -23,6 +23,8 @@
#include "domain.h"
#include "style_region.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "force.h"
#include "kspace.h"
#include "update.h"
@ -551,7 +553,8 @@ void Domain::pbc()
void Domain::image_check()
{
int i,j,k;
int i,j,k,n,imol,iatom;
tagint tagprev;
// only need to check if system is molecular and some dimension is periodic
// if running verlet/split, don't check on KSpace partition since
@ -581,8 +584,15 @@ void Domain::image_check()
// flag if any bond component is longer than non-periodic box length
// which means image flags in that dimension were different
int molecular = atom->molecular;
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
double delx,dely,delz;
@ -590,9 +600,25 @@ void Domain::image_check()
int nmissing = 0;
int flag = 0;
for (i = 0; i < nlocal; i++)
for (j = 0; j < num_bond[i]; j++) {
for (i = 0; i < nlocal; i++) {
if (molecular == 1) n = num_bond[i];
else {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
n = onemols[imol]->num_bond[iatom];
}
for (j = 0; j < n; j++) {
if (molecular == 1) {
if (bond_type[i][j] <= 0) continue;
k = atom->map(bond_atom[i][j]);
} else {
if (onemols[imol]->bond_type[iatom][j] < 0) continue;
tagprev = tag[i] - iatom - 1;
k = atom->map(onemols[imol]->bond_atom[iatom][j]+tagprev);
}
if (k == -1) {
nmissing++;
if (lostbond == ERROR)
@ -611,6 +637,7 @@ void Domain::image_check()
if (!yperiodic && dely > yprd) flag = 1;
if (dimension == 3 && !zperiodic && delz > zprd) flag = 1;
}
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
@ -636,7 +663,8 @@ void Domain::image_check()
void Domain::box_too_small_check()
{
int i,j,k;
int i,j,k,n,imol,iatom;
tagint tagprev;
// only need to check if system is molecular and some dimension is periodic
// if running verlet/split, don't check on KSpace partition since
@ -653,10 +681,16 @@ void Domain::box_too_small_check()
// in this case, image_check() should warn,
// assuming 2 atoms have consistent image flags
int molecular = atom->molecular;
double **x = atom->x;
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
double **x = atom->x;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int nlocal = atom->nlocal;
double delx,dely,delz,rsq,r;
@ -665,16 +699,32 @@ void Domain::box_too_small_check()
int lostbond = output->thermo->lostbond;
int nmissing = 0;
for (i = 0; i < nlocal; i++)
for (j = 0; j < num_bond[i]; j++) {
for (i = 0; i < nlocal; i++) {
if (molecular == 1) n = num_bond[i];
else {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
n = onemols[imol]->num_bond[iatom];
}
for (j = 0; j < n; j++) {
if (molecular == 1) {
if (bond_type[i][j] <= 0) continue;
k = atom->map(bond_atom[i][j]);
} else {
if (onemols[imol]->bond_type[iatom][j] < 0) continue;
tagprev = tag[i] - iatom - 1;
k = atom->map(onemols[imol]->bond_atom[iatom][j]+tagprev);
}
if (k == -1) {
nmissing++;
if (lostbond == ERROR)
error->one(FLERR,"Bond atom missing in box size check");
continue;
}
delx = x[i][0] - x[k][0];
dely = x[i][1] - x[k][1];
delz = x[i][2] - x[k][2];
@ -682,6 +732,7 @@ void Domain::box_too_small_check()
rsq = delx*delx + dely*dely + delz*delz;
maxbondme = MAX(maxbondme,rsq);
}
}
if (lostbond == WARN) {
int all;

4
src/domain.h
View File

@ -131,8 +131,10 @@ class Domain : protected Pointers {
// minimum image convention check
// return 1 if any distance > 1/2 of box size
// indicates a special neighbor is actually not in a bond,
// but is a far-away image that should be treated as an unbonded neighbor
// inline since called from neighbor build inner loop
//
inline int minimum_image_check(double dx, double dy, double dz) {
if (xperiodic && fabs(dx) > xprd_half) return 1;
if (yperiodic && fabs(dy) > yprd_half) return 1;

3
src/dump.cpp
View File

@ -106,6 +106,9 @@ Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
char *ptr;
if (ptr = strchr(filename,'%')) {
if (strstr(style,"mpiio"))
error->all(FLERR,
"Dump file MPI-IO output not allowed with '%' in filename");
multiproc = 1;
nclusterprocs = 1;
filewriter = 1;

33
src/dump_image.cpp
View File

@ -18,6 +18,8 @@
#include "dump_image.h"
#include "image.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "group.h"
#include "force.h"
@ -646,7 +648,8 @@ void DumpImage::view_params()
void DumpImage::create_image()
{
int i,j,m,itype,atom1,atom2;
int i,j,m,n,itype,atom1,atom2,imol,iatom,btype;
tagint tagprev;
double diameter,delx,dely,delz;
double *color,*color1,*color2;
double xmid[3];
@ -700,10 +703,14 @@ void DumpImage::create_image()
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
int *num_bond = atom->num_bond;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
int newton_bond = force->newton_bond;
int molecular = atom->molecular;
Molecule **onemols = atom->avec->onemols;
// communicate choose flag for ghost atoms to know if they are selected
// if bcolor/bdiam = ATOM, setup bufcopy to comm atom color/diam attributes
@ -735,11 +742,27 @@ void DumpImage::create_image()
for (i = 0; i < nchoose; i++) {
atom1 = clist[i];
for (m = 0; m < num_bond[atom1]; m++) {
if (molecular == 1) n = num_bond[atom1];
else {
if (molindex[atom1] < 0) continue;
imol = molindex[atom1];
iatom = molatom[atom1];
n = onemols[imol]->num_bond[iatom];
}
for (m = 0; m < n; m++) {
if (molecular == 1) {
btype = bond_type[atom1][m];
atom2 = atom->map(bond_atom[atom1][m]);
} else {
tagprev = tag[i] - iatom - 1;
btype = atom->map(onemols[imol]->bond_type[atom1][m]);
atom2 = atom->map(onemols[imol]->bond_atom[iatom][m]+tagprev);
}
if (atom2 < 0 || !chooseghost[atom2]) continue;
if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
if (bond_type[atom1][m] == 0) continue;
if (btype == 0) continue;
if (bcolor == ATOM) {
if (acolor == TYPE) {
@ -753,7 +776,7 @@ void DumpImage::create_image()
color2 = image->map_value2color(0,bufcopy[atom2][0]);
}
} else if (bcolor == TYPE) {
itype = bond_type[atom1][m];
itype = btype;
if (itype < 0) itype = -itype;
color = bcolortype[itype];
}
@ -771,7 +794,7 @@ void DumpImage::create_image()
diameter = MIN(bufcopy[atom1][1],bufcopy[atom2][1]);
}
} else if (bdiam == TYPE) {
itype = bond_type[atom1][m];
itype = btype;
if (itype < 0) itype = -itype;
diameter = bdiamtype[itype];
}

24
src/finish.cpp
View File

@ -19,6 +19,8 @@
#include "timer.h"
#include "universe.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "comm.h"
#include "force.h"
#include "kspace.h"
@ -608,10 +610,26 @@ void Finish::end(int flag)
int nspec;
double nspec_all;
if (atom->molecular) {
nspec = 0;
if (atom->molecular == 1) {
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) nspec += atom->nspecial[i][2];
nspec = 0;
for (i = 0; i < nlocal; i++) nspec += nspecial[i][2];
tmp = nspec;
MPI_Allreduce(&tmp,&nspec_all,1,MPI_DOUBLE,MPI_SUM,world);
} else if (atom->molecular == 2) {
Molecule **onemols = atom->avec->onemols;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int nlocal = atom->nlocal;
int imol,iatom;
nspec = 0;
for (i = 0; i < nlocal; i++) {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
nspec += onemols[imol]->nspecial[iatom][2];
}
tmp = nspec;
MPI_Allreduce(&tmp,&nspec_all,1,MPI_DOUBLE,MPI_SUM,world);
}

2
src/input.cpp
View File

@ -1504,7 +1504,7 @@ void Input::special_bonds()
// if simulation box defined and saved values changed, redo special list
if (domain->box_exist && atom->molecular) {
if (domain->box_exist && atom->molecular == 1) {
if (lj2 != force->special_lj[2] || lj3 != force->special_lj[3] ||
coul2 != force->special_coul[2] || coul3 != force->special_coul[3] ||
angle != force->special_angle ||

152
src/molecule.cpp
View File

@ -348,17 +348,6 @@ void Molecule::read(int flag)
if (flag == 0) {
if (natoms == 0) error->all(FLERR,"No atom count in molecule file");
if (nbonds && !atom->avec->bonds_allow)
error->all(FLERR,"Bonds in molecule file not supported by atom style");
if (nangles && !atom->avec->angles_allow)
error->all(FLERR,"Angles in molecule file not supported by atom style");
if (ndihedrals && !atom->avec->dihedrals_allow)
error->all(FLERR,
"Dihedrals in molecule file not supported by atom style");
if (nimpropers && !atom->avec->impropers_allow)
error->all(FLERR,
"Impropers in molecule file not supported by atom style");
}
// count = vector for tallying bonds,angles,etc per atom
@ -455,32 +444,6 @@ void Molecule::read(int flag)
// error check
if (flag == 0) {
if (qflag && !atom->q_flag)
error->all(FLERR,"Molecule file has undefined atom property");
if (radiusflag && !atom->radius_flag)
error->all(FLERR,"Molecule file has undefined atom property");
if (rmassflag && !atom->rmass_flag)
error->all(FLERR,"Molecule file has undefined atom property");
if (bondflag && !atom->avec->bonds_allow)
error->all(FLERR,"Invalid molecule file section: Bonds");
if (angleflag && !atom->avec->angles_allow)
error->all(FLERR,"Invalid molecule file section: Angles");
if (dihedralflag && !atom->avec->dihedrals_allow)
error->all(FLERR,"Invalid molecule file section: Dihedrals");
if (improperflag && !atom->avec->impropers_allow)
error->all(FLERR,"Invalid molecule file section: Impropers");
if (bond_per_atom > atom->bond_per_atom ||
angle_per_atom > atom->angle_per_atom ||
dihedral_per_atom > atom->dihedral_per_atom ||
improper_per_atom > atom->improper_per_atom)
error->all(FLERR,"Molecule file bond/angle/etc counts "
"per atom are too large");
// test for maxspecial > atom->maxspecial is done when molecules added
// in Atom::add_molecule_atom()
if ((nspecialflag && !specialflag) || (!nspecialflag && specialflag))
error->all(FLERR,"Molecule file needs both Special Bond sections");
if (specialflag && !bondflag)
@ -512,7 +475,7 @@ void Molecule::coords(char *line)
/* ----------------------------------------------------------------------
read types from file
set maxtype = max of any atom type
set ntypes = max of any atom type
------------------------------------------------------------------------- */
void Molecule::types(char *line)
@ -524,11 +487,11 @@ void Molecule::types(char *line)
}
for (int i = 0; i < natoms; i++)
if (type[i] <= 0 || type[i] > atom->ntypes)
if (type[i] <= 0)
error->all(FLERR,"Invalid atom type in molecule file");
for (int i = 0; i < natoms; i++)
maxtype = MAX(maxtype,type[i]);
ntypes = MAX(ntypes,type[i]);
}
/* ----------------------------------------------------------------------
@ -580,6 +543,7 @@ void Molecule::masses(char *line)
/* ----------------------------------------------------------------------
read bonds from file
set nbondtypes = max type of any bond
store each with both atoms if newton_bond = 0
if flag = 0, just count bonds/atom
if flag = 1, store them with atoms
@ -602,11 +566,12 @@ void Molecule::bonds(int flag, char *line)
if (atom1 <= 0 || atom1 > natoms ||
atom2 <= 0 || atom2 > natoms)
error->one(FLERR,"Invalid atom ID in Bonds section of molecule file");
if (itype <= 0 || itype > atom->nbondtypes)
if (itype <= 0)
error->one(FLERR,"Invalid bond type in Bonds section of molecule file");
if (flag) {
m = atom1-1;
nbondtypes = MAX(nbondtypes,itype);
bond_type[m][num_bond[m]] = itype;
bond_atom[m][num_bond[m]] = atom2;
num_bond[m]++;
@ -656,11 +621,12 @@ void Molecule::angles(int flag, char *line)
atom2 <= 0 || atom2 > natoms ||
atom3 <= 0 || atom3 > natoms)
error->one(FLERR,"Invalid atom ID in Angles section of molecule file");
if (itype <= 0 || itype > atom->nangletypes)
if (itype <= 0)
error->one(FLERR,"Invalid angle type in Angles section of molecule file");
if (flag) {
m = atom2-1;
nangletypes = MAX(nangletypes,itype);
angle_type[m][num_angle[m]] = itype;
angle_atom1[m][num_angle[m]] = atom1;
angle_atom2[m][num_angle[m]] = atom2;
@ -725,12 +691,13 @@ void Molecule::dihedrals(int flag, char *line)
atom4 <= 0 || atom4 > natoms)
error->one(FLERR,
"Invalid atom ID in dihedrals section of molecule file");
if (itype <= 0 || itype > atom->ndihedraltypes)
if (itype <= 0)
error->one(FLERR,
"Invalid dihedral type in dihedrals section of molecule file");
if (flag) {
m = atom2-1;
ndihedraltypes = MAX(ndihedraltypes,itype);
dihedral_type[m][num_dihedral[m]] = itype;
dihedral_atom1[m][num_dihedral[m]] = atom1;
dihedral_atom2[m][num_dihedral[m]] = atom2;
@ -802,12 +769,13 @@ void Molecule::impropers(int flag, char *line)
atom4 <= 0 || atom4 > natoms)
error->one(FLERR,
"Invalid atom ID in impropers section of molecule file");
if (itype <= 0 || itype > atom->nimpropertypes)
if (itype <= 0)
error->one(FLERR,
"Invalid improper type in impropers section of molecule file");
if (flag) {
m = atom2-1;
nimpropertypes = MAX(nimpropertypes,itype);
improper_type[m][num_improper[m]] = itype;
improper_atom1[m][num_improper[m]] = atom1;
improper_atom2[m][num_improper[m]] = atom2;
@ -977,14 +945,79 @@ void Molecule::shaketype_read(char *line)
int m = shake_flag[i];
if (m == 1) m = 3;
for (int j = 0; j < m-1; j++)
if (shake_type[i][j] <= 0 || shake_type[i][j] > atom->nbondtypes)
if (shake_type[i][j] <= 0)
error->all(FLERR,"Invalid shake bond type in molecule file");
if (shake_flag[i] == 1)
if (shake_type[i][2] <= 0 || shake_type[i][2] > atom->nangletypes)
if (shake_type[i][2] <= 0)
error->all(FLERR,"Invalid shake angle type in molecule file");
}
}
/* ----------------------------------------------------------------------
error check molecule attributes and topology against system settings
flag = 0, just check this molecule
flag = 1, check all molecules in set, this is 1st molecule in set
------------------------------------------------------------------------- */
void Molecule::check_attributes(int flag)
{
int n = 1;
if (flag) n = nset;
int imol = atom->find_molecule(id);
for (int i = imol; i < imol+n; i++) {
Molecule *onemol = atom->molecules[imol];
// check per-atom attributes of molecule
// warn if not a match
int mismatch = 0;
if (onemol->qflag && !atom->q_flag) mismatch = 1;
if (onemol->radiusflag && !atom->radius_flag) mismatch = 1;
if (onemol->rmassflag && !atom->rmass_flag) mismatch = 1;
if (mismatch && me == 0)
error->warning(FLERR,
"Molecule attributes do not match system attributes");
// for all atom styles, check nbondtype,etc
mismatch = 0;
if (atom->nbondtypes < onemol->nbondtypes) mismatch = 1;
if (atom->nangletypes < onemol->nangletypes) mismatch = 1;
if (atom->ndihedraltypes < onemol->ndihedraltypes) mismatch = 1;
if (atom->nimpropertypes < onemol->nimpropertypes) mismatch = 1;
if (mismatch)
error->all(FLERR,"Molecule topology type exceeds system topology type");
// for molecular atom styles, check bond_per_atom,etc + maxspecial
// do not check for atom style template, since nothing stored per atom
if (atom->molecular == 1) {
if (atom->avec->bonds_allow &&
atom->bond_per_atom < onemol->bond_per_atom) mismatch = 1;
if (atom->avec->angles_allow &&
atom->angle_per_atom < onemol->angle_per_atom) mismatch = 1;
if (atom->avec->dihedrals_allow &&
atom->dihedral_per_atom < onemol->dihedral_per_atom) mismatch = 1;
if (atom->avec->impropers_allow &&
atom->improper_per_atom < onemol->improper_per_atom) mismatch = 1;
if (atom->maxspecial < onemol->maxspecial) mismatch = 1;
if (mismatch)
error->all(FLERR,"Molecule toplogy/atom exceeds system topology/atom");
}
// warn if molecule topology defined but no special settings
if (onemol->bondflag && !onemol->specialflag)
if (me == 0) error->warning(FLERR,"Molecule has bond topology "
"but no special bond settings");
}
}
/* ----------------------------------------------------------------------
init all data structures to empty
------------------------------------------------------------------------- */
@ -993,7 +1026,9 @@ void Molecule::initialize()
{
natoms = 0;
nbonds = nangles = ndihedrals = nimpropers = 0;
maxtype = 0;
ntypes = 0;
nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
maxspecial = 0;
@ -1040,6 +1075,7 @@ void Molecule::initialize()
/* ----------------------------------------------------------------------
allocate all data structures
also initialize values for data structures that are always allocated
------------------------------------------------------------------------- */
void Molecule::allocate()
@ -1050,8 +1086,23 @@ void Molecule::allocate()
if (radiusflag) memory->create(radius,natoms,"molecule:radius");
if (rmassflag) memory->create(rmass,natoms,"molecule:rmass");
if (bondflag) {
// always allocate num_bond,num_angle,etc and nspecial even if not in file
// initialize to 0 even if not in molecule file
// this is so methods that use these arrays don't have to check they exist
memory->create(num_bond,natoms,"molecule:num_bond");
for (int i = 0; i < natoms; i++) num_bond[i] = 0;
memory->create(num_angle,natoms,"molecule:num_angle");
for (int i = 0; i < natoms; i++) num_angle[i] = 0;
memory->create(num_dihedral,natoms,"molecule:num_dihedral");
for (int i = 0; i < natoms; i++) num_dihedral[i] = 0;
memory->create(num_improper,natoms,"molecule:num_improper");
for (int i = 0; i < natoms; i++) num_improper[i] = 0;
memory->create(nspecial,natoms,3,"molecule:nspecial");
for (int i = 0; i < natoms; i++)
nspecial[i][0] = nspecial[i][1] = nspecial[i][2] = 0;
if (bondflag) {
memory->create(bond_type,natoms,bond_per_atom,
"molecule:bond_type");
memory->create(bond_atom,natoms,bond_per_atom,
@ -1059,7 +1110,6 @@ void Molecule::allocate()
}
if (angleflag) {
memory->create(num_angle,natoms,"molecule:num_angle");
memory->create(angle_type,natoms,angle_per_atom,
"molecule:angle_type");
memory->create(angle_atom1,natoms,angle_per_atom,
@ -1071,7 +1121,6 @@ void Molecule::allocate()
}
if (dihedralflag) {
memory->create(num_dihedral,natoms,"molecule:num_dihedral");
memory->create(dihedral_type,natoms,dihedral_per_atom,
"molecule:dihedral_type");
memory->create(dihedral_atom1,natoms,dihedral_per_atom,
@ -1085,7 +1134,6 @@ void Molecule::allocate()
}
if (improperflag) {
memory->create(num_improper,natoms,"molecule:num_improper");
memory->create(improper_type,natoms,improper_per_atom,
"molecule:improper_type");
memory->create(improper_atom1,natoms,improper_per_atom,
@ -1098,8 +1146,6 @@ void Molecule::allocate()
"molecule:improper_atom4");
}
if (nspecialflag)
memory->create(nspecial,natoms,3,"molecule:nspecial");
if (specialflag)
memory->create(special,natoms,maxspecial,"molecule:special");

4
src/molecule.h
View File

@ -28,10 +28,11 @@ class Molecule : protected Pointers {
int natoms;
int nbonds,nangles,ndihedrals,nimpropers;
int ntypes;
int nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
// max bond,angle,etc per atom
int maxtype;
int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
int maxspecial;
@ -105,6 +106,7 @@ class Molecule : protected Pointers {
void compute_mass();
void compute_com();
void compute_inertia();
void check_attributes(int);
private:
int me;

336
src/neigh_bond.cpp
View File

@ -13,6 +13,8 @@
#include "neighbor.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "force.h"
#include "update.h"
#include "domain.h"
@ -92,6 +94,78 @@ void Neighbor::bond_all()
/* ---------------------------------------------------------------------- */
void Neighbor::bond_template()
{
int i,m,atom1;
int imol,iatom;
tagint tagprev;
int *num_bond;
int **bond_atom,**bond_type;
Molecule **onemols = atom->avec->onemols;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
int lostbond = output->thermo->lostbond;
int nmissing = 0;
nbondlist = 0;
for (i = 0; i < nlocal; i++) {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
num_bond = onemols[imol]->num_bond;
bond_atom = onemols[imol]->bond_atom;
bond_type = onemols[imol]->bond_type;
for (m = 0; m < num_bond[iatom]; m++) {
if (bond_type[iatom][m] <= 0) continue;
atom1 = atom->map(bond_atom[iatom][m]+tagprev);
if (atom1 == -1) {
nmissing++;
if (lostbond == ERROR) {
char str[128];
sprintf(str,"Bond atoms " TAGINT_FORMAT " " TAGINT_FORMAT
" missing on proc %d at step " BIGINT_FORMAT,
tag[i],bond_atom[iatom][m]+tagprev,me,update->ntimestep);
error->one(FLERR,str);
}
continue;
}
atom1 = domain->closest_image(i,atom1);
if (newton_bond || i < atom1) {
if (nbondlist == maxbond) {
maxbond += BONDDELTA;
memory->grow(bondlist,maxbond,3,"neighbor:bondlist");
}
bondlist[nbondlist][0] = i;
bondlist[nbondlist][1] = atom1;
bondlist[nbondlist][2] = bond_type[iatom][m];
nbondlist++;
}
}
}
if (cluster_check) bond_check();
if (lostbond == IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Bond atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
}
/* ---------------------------------------------------------------------- */
void Neighbor::bond_partial()
{
int i,m,atom1;
@ -240,6 +314,88 @@ void Neighbor::angle_all()
/* ---------------------------------------------------------------------- */
void Neighbor::angle_template()
{
int i,m,atom1,atom2,atom3;
int imol,iatom;
tagint tagprev;
int *num_angle;
int **angle_atom1,**angle_atom2,**angle_atom3,**angle_type;
Molecule **onemols = atom->avec->onemols;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
int lostbond = output->thermo->lostbond;
int nmissing = 0;
nanglelist = 0;
for (i = 0; i < nlocal; i++) {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
num_angle = onemols[imol]->num_angle;
angle_atom1 = onemols[imol]->angle_atom1;
angle_atom2 = onemols[imol]->angle_atom2;
angle_atom3 = onemols[imol]->angle_atom3;
angle_type = onemols[imol]->angle_type;
for (m = 0; m < num_angle[iatom]; m++) {
if (angle_type[iatom][m] <= 0) continue;
atom1 = atom->map(angle_atom1[iatom][m]+tagprev);
atom2 = atom->map(angle_atom2[iatom][m]+tagprev);
atom3 = atom->map(angle_atom3[iatom][m]+tagprev);
if (atom1 == -1 || atom2 == -1 || atom3 == -1) {
nmissing++;
if (lostbond == ERROR) {
char str[128];
sprintf(str,"Angle atoms "
TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
" missing on proc %d at step " BIGINT_FORMAT,
angle_atom1[iatom][m]+tagprev,angle_atom2[iatom][m]+tagprev,
angle_atom3[iatom][m]+tagprev,
me,update->ntimestep);
error->one(FLERR,str);
}
continue;
}
atom1 = domain->closest_image(i,atom1);
atom2 = domain->closest_image(i,atom2);
atom3 = domain->closest_image(i,atom3);
if (newton_bond || (i <= atom1 && i <= atom2 && i <= atom3)) {
if (nanglelist == maxangle) {
maxangle += BONDDELTA;
memory->grow(anglelist,maxangle,4,"neighbor:anglelist");
}
anglelist[nanglelist][0] = atom1;
anglelist[nanglelist][1] = atom2;
anglelist[nanglelist][2] = atom3;
anglelist[nanglelist][3] = angle_type[iatom][m];
nanglelist++;
}
}
}
if (cluster_check) angle_check();
if (lostbond == IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Angle atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
}
/* ---------------------------------------------------------------------- */
void Neighbor::angle_partial()
{
int i,m,atom1,atom2,atom3;
@ -417,6 +573,96 @@ void Neighbor::dihedral_all()
/* ---------------------------------------------------------------------- */
void Neighbor::dihedral_template()
{
int i,m,atom1,atom2,atom3,atom4;
int imol,iatom;
tagint tagprev;
int *num_dihedral;
int **dihedral_atom1,**dihedral_atom2,**dihedral_atom3,**dihedral_atom4;
int **dihedral_type;
Molecule **onemols = atom->avec->onemols;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
int lostbond = output->thermo->lostbond;
int nmissing = 0;
ndihedrallist = 0;
for (i = 0; i < nlocal; i++) {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
num_dihedral = onemols[imol]->num_dihedral;
dihedral_atom1 = onemols[imol]->dihedral_atom1;
dihedral_atom2 = onemols[imol]->dihedral_atom2;
dihedral_atom3 = onemols[imol]->dihedral_atom3;
dihedral_atom4 = onemols[imol]->dihedral_atom4;
dihedral_type = onemols[imol]->dihedral_type;
for (m = 0; m < num_dihedral[iatom]; m++) {
atom1 = atom->map(dihedral_atom1[iatom][m]+tagprev);
atom2 = atom->map(dihedral_atom2[iatom][m]+tagprev);
atom3 = atom->map(dihedral_atom3[iatom][m]+tagprev);
atom4 = atom->map(dihedral_atom4[iatom][m]+tagprev);
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1) {
nmissing++;
if (lostbond == ERROR) {
char str[128];
sprintf(str,"Dihedral atoms "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT
" missing on proc %d at step " BIGINT_FORMAT,
dihedral_atom1[iatom][m]+tagprev,
dihedral_atom2[iatom][m]+tagprev,
dihedral_atom3[iatom][m]+tagprev,
dihedral_atom4[iatom][m]+tagprev,
me,update->ntimestep);
error->one(FLERR,str);
}
continue;
}
atom1 = domain->closest_image(i,atom1);
atom2 = domain->closest_image(i,atom2);
atom3 = domain->closest_image(i,atom3);
atom4 = domain->closest_image(i,atom4);
if (newton_bond ||
(i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)) {
if (ndihedrallist == maxdihedral) {
maxdihedral += BONDDELTA;
memory->grow(dihedrallist,maxdihedral,5,"neighbor:dihedrallist");
}
dihedrallist[ndihedrallist][0] = atom1;
dihedrallist[ndihedrallist][1] = atom2;
dihedrallist[ndihedrallist][2] = atom3;
dihedrallist[ndihedrallist][3] = atom4;
dihedrallist[ndihedrallist][4] = dihedral_type[iatom][m];
ndihedrallist++;
}
}
}
if (cluster_check) dihedral_check(ndihedrallist,dihedrallist);
if (lostbond == IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Dihedral atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
}
/* ---------------------------------------------------------------------- */
void Neighbor::dihedral_partial()
{
int i,m,atom1,atom2,atom3,atom4;
@ -618,6 +864,96 @@ void Neighbor::improper_all()
/* ---------------------------------------------------------------------- */
void Neighbor::improper_template()
{
int i,m,atom1,atom2,atom3,atom4;
int imol,iatom;
tagint tagprev;
int *num_improper;
int **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4;
int **improper_type;
Molecule **onemols = atom->avec->onemols;
tagint *tag = atom->tag;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
int lostbond = output->thermo->lostbond;
int nmissing = 0;
nimproperlist = 0;
for (i = 0; i < nlocal; i++) {
if (molindex[i] < 0) continue;
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
num_improper = onemols[imol]->num_improper;
improper_atom1 = onemols[imol]->improper_atom1;
improper_atom2 = onemols[imol]->improper_atom2;
improper_atom3 = onemols[imol]->improper_atom3;
improper_atom4 = onemols[imol]->improper_atom4;
improper_type = onemols[imol]->improper_type;
for (m = 0; m < num_improper[iatom]; m++) {
atom1 = atom->map(improper_atom1[iatom][m]+tagprev);
atom2 = atom->map(improper_atom2[iatom][m]+tagprev);
atom3 = atom->map(improper_atom3[iatom][m]+tagprev);
atom4 = atom->map(improper_atom4[iatom][m]+tagprev);
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1) {
nmissing++;
if (lostbond == ERROR) {
char str[128];
sprintf(str,"Improper atoms "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT
" missing on proc %d at step " BIGINT_FORMAT,
improper_atom1[iatom][m]+tagprev,
improper_atom2[iatom][m]+tagprev,
improper_atom3[iatom][m]+tagprev,
improper_atom4[iatom][m]+tagprev,
me,update->ntimestep);
error->one(FLERR,str);
}
continue;
}
atom1 = domain->closest_image(i,atom1);
atom2 = domain->closest_image(i,atom2);
atom3 = domain->closest_image(i,atom3);
atom4 = domain-> closest_image(i,atom4);
if (newton_bond ||
(i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)) {
if (nimproperlist == maximproper) {
maximproper += BONDDELTA;
memory->grow(improperlist,maximproper,5,"neighbor:improperlist");
}
improperlist[nimproperlist][0] = atom1;
improperlist[nimproperlist][1] = atom2;
improperlist[nimproperlist][2] = atom3;
improperlist[nimproperlist][3] = atom4;
improperlist[nimproperlist][4] = improper_type[iatom][m];
nimproperlist++;
}
}
}
if (cluster_check) dihedral_check(nimproperlist,improperlist);
if (lostbond == IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Improper atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
}
/* ---------------------------------------------------------------------- */
void Neighbor::improper_partial()
{
int i,m,atom1,atom2,atom3,atom4;

157
src/neigh_full.cpp
View File

@ -14,6 +14,8 @@
#include "neighbor.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "group.h"
@ -28,26 +30,32 @@ using namespace LAMMPS_NS;
void Neighbor::full_nsq(NeighList *list)
{
int i,j,n,itype,jtype,which,bitmask;
int i,j,n,itype,jtype,which,bitmask,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int molecular = atom->molecular;
if (includegroup) {
nlocal = atom->nfirst;
bitmask = group->bitmask[includegroup];
}
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -66,6 +74,11 @@ void Neighbor::full_nsq(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms, owned and ghost
// skip i = j
@ -82,7 +95,13 @@ void Neighbor::full_nsq(NeighList *list)
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -111,21 +130,27 @@ void Neighbor::full_nsq(NeighList *list)
void Neighbor::full_nsq_ghost(NeighList *list)
{
int i,j,n,itype,jtype,which;
int i,j,n,itype,jtype,which,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
@ -145,6 +170,11 @@ void Neighbor::full_nsq_ghost(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms, owned and ghost
// skip i = j
@ -162,7 +192,13 @@ void Neighbor::full_nsq_ghost(NeighList *list)
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -204,7 +240,8 @@ void Neighbor::full_nsq_ghost(NeighList *list)
void Neighbor::full_bin(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which;
int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
@ -212,18 +249,23 @@ void Neighbor::full_bin(NeighList *list)
bin_atoms();
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -244,6 +286,11 @@ void Neighbor::full_bin(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in surrounding bins in stencil including self
// skip i = j
@ -264,7 +311,13 @@ void Neighbor::full_bin(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -294,7 +347,8 @@ void Neighbor::full_bin(NeighList *list)
void Neighbor::full_bin_ghost(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which;
int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int xbin,ybin,zbin,xbin2,ybin2,zbin2;
int *neighptr;
@ -303,17 +357,22 @@ void Neighbor::full_bin_ghost(NeighList *list)
bin_atoms();
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
@ -336,6 +395,11 @@ void Neighbor::full_bin_ghost(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in surrounding bins in stencil including self
// when i is a ghost atom, must check if stencil bin is out of bounds
@ -358,7 +422,13 @@ void Neighbor::full_bin_ghost(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -413,7 +483,8 @@ void Neighbor::full_bin_ghost(NeighList *list)
void Neighbor::full_multi(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which,ns;
int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
double *cutsq,*distsq;
@ -422,20 +493,23 @@ void Neighbor::full_multi(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -447,6 +521,8 @@ void Neighbor::full_multi(NeighList *list)
int inum = 0;
ipage->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
@ -455,6 +531,11 @@ void Neighbor::full_multi(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in other bins in stencil, including self
// skip if i,j neighbor cutoff is less than bin distance
@ -480,7 +561,13 @@ void Neighbor::full_multi(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;

142
src/neigh_half_bin.cpp
View File

@ -14,6 +14,8 @@
#include "neighbor.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
@ -29,7 +31,8 @@ using namespace LAMMPS_NS;
void Neighbor::half_bin_no_newton(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which;
int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
@ -37,19 +40,22 @@ void Neighbor::half_bin_no_newton(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
@ -61,6 +67,8 @@ void Neighbor::half_bin_no_newton(NeighList *list)
int inum = 0;
ipage->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
@ -69,6 +77,11 @@ void Neighbor::half_bin_no_newton(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in other bins in stencil including self
// only store pair if i < j
@ -91,11 +104,18 @@ void Neighbor::half_bin_no_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
// OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
@ -123,7 +143,8 @@ void Neighbor::half_bin_no_newton(NeighList *list)
void Neighbor::half_bin_no_newton_ghost(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which;
int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int xbin,ybin,zbin,xbin2,ybin2,zbin2;
int *neighptr;
@ -132,19 +153,23 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
@ -157,6 +182,8 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list)
int inum = 0;
ipage->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nall; i++) {
n = 0;
neighptr = ipage->vget();
@ -165,6 +192,11 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in other bins in stencil including self
// when i is a ghost atom, must check if stencil bin is out of bounds
@ -191,7 +223,13 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -246,7 +284,8 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list)
void Neighbor::half_bin_newton(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which;
int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
@ -254,19 +293,22 @@ void Neighbor::half_bin_newton(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
@ -278,6 +320,8 @@ void Neighbor::half_bin_newton(NeighList *list)
int inum = 0;
ipage->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
@ -286,6 +330,11 @@ void Neighbor::half_bin_newton(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// 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
@ -310,7 +359,13 @@ void Neighbor::half_bin_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -335,7 +390,13 @@ void Neighbor::half_bin_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -365,7 +426,8 @@ void Neighbor::half_bin_newton(NeighList *list)
void Neighbor::half_bin_newton_tri(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which;
int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
@ -373,19 +435,22 @@ void Neighbor::half_bin_newton_tri(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
@ -397,6 +462,8 @@ void Neighbor::half_bin_newton_tri(NeighList *list)
int inum = 0;
ipage->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
@ -405,6 +472,11 @@ void Neighbor::half_bin_newton_tri(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
@ -434,7 +506,13 @@ void Neighbor::half_bin_newton_tri(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;

113
src/neigh_half_multi.cpp
View File

@ -14,6 +14,8 @@
#include "neighbor.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
@ -30,7 +32,8 @@ using namespace LAMMPS_NS;
void Neighbor::half_multi_no_newton(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which,ns;
int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
double *cutsq,*distsq;
@ -39,20 +42,23 @@ void Neighbor::half_multi_no_newton(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -64,6 +70,8 @@ void Neighbor::half_multi_no_newton(NeighList *list)
int inum = 0;
ipage->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
@ -72,6 +80,11 @@ void Neighbor::half_multi_no_newton(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in other bins in stencil including self
// only store pair if i < j
@ -99,7 +112,13 @@ void Neighbor::half_multi_no_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -129,7 +148,8 @@ void Neighbor::half_multi_no_newton(NeighList *list)
void Neighbor::half_multi_newton(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which,ns;
int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
double *cutsq,*distsq;
@ -138,20 +158,23 @@ void Neighbor::half_multi_newton(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -163,6 +186,8 @@ void Neighbor::half_multi_newton(NeighList *list)
int inum = 0;
ipage->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
@ -171,6 +196,11 @@ void Neighbor::half_multi_newton(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// 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
@ -195,7 +225,13 @@ void Neighbor::half_multi_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -226,7 +262,13 @@ void Neighbor::half_multi_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -256,7 +298,8 @@ void Neighbor::half_multi_newton(NeighList *list)
void Neighbor::half_multi_newton_tri(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which,ns;
int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
double *cutsq,*distsq;
@ -265,20 +308,23 @@ void Neighbor::half_multi_newton_tri(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -290,6 +336,8 @@ void Neighbor::half_multi_newton_tri(NeighList *list)
int inum = 0;
ipage->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
@ -298,6 +346,11 @@ void Neighbor::half_multi_newton_tri(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in bins, including self, in stencil
// skip if i,j neighbor cutoff is less than bin distance
@ -334,7 +387,13 @@ void Neighbor::half_multi_newton_tri(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;

99
src/neigh_half_nsq.cpp
View File

@ -14,6 +14,8 @@
#include "neighbor.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "group.h"
#include "error.h"
@ -28,26 +30,32 @@ using namespace LAMMPS_NS;
void Neighbor::half_nsq_no_newton(NeighList *list)
{
int i,j,n,itype,jtype,which,bitmask;
int i,j,n,itype,jtype,which,bitmask,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int molecular = atom->molecular;
if (includegroup) {
nlocal = atom->nfirst;
bitmask = group->bitmask[includegroup];
}
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -66,6 +74,11 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over remaining atoms, owned and ghost
// only store pair if i < j
@ -82,7 +95,13 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -112,26 +131,32 @@ void Neighbor::half_nsq_no_newton(NeighList *list)
void Neighbor::half_nsq_no_newton_ghost(NeighList *list)
{
int i,j,n,itype,jtype,which,bitmask;
int i,j,n,itype,jtype,which,bitmask,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int molecular = atom->molecular;
if (includegroup) {
nlocal = atom->nfirst;
bitmask = group->bitmask[includegroup];
}
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -150,6 +175,11 @@ void Neighbor::half_nsq_no_newton_ghost(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over remaining atoms, owned and ghost
// only store pair if i < j
@ -171,7 +201,13 @@ void Neighbor::half_nsq_no_newton_ghost(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -215,28 +251,32 @@ void Neighbor::half_nsq_no_newton_ghost(NeighList *list)
void Neighbor::half_nsq_newton(NeighList *list)
{
int i,j,n,itype,jtype,itag,jtag,which,bitmask;
int i,j,n,itype,jtype,itag,jtag,which,bitmask,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int molecular = atom->molecular;
if (includegroup) {
nlocal = atom->nfirst;
bitmask = group->bitmask[includegroup];
}
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -245,6 +285,8 @@ void Neighbor::half_nsq_newton(NeighList *list)
int inum = 0;
ipage->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
@ -254,6 +296,11 @@ void Neighbor::half_nsq_newton(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over remaining atoms, owned and ghost
// itag = jtag is possible for long cutoffs that include images of self
@ -286,7 +333,13 @@ void Neighbor::half_nsq_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;

182
src/neigh_respa.cpp
View File

@ -14,6 +14,8 @@
#include "neighbor.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "group.h"
#include "my_page.h"
@ -30,28 +32,32 @@ using namespace LAMMPS_NS;
void Neighbor::respa_nsq_no_newton(NeighList *list)
{
int i,j,n,itype,jtype,n_inner,n_middle,bitmask;
int i,j,n,itype,jtype,n_inner,n_middle,bitmask,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int molecular = atom->molecular;
if (includegroup) {
nlocal = atom->nfirst;
bitmask = group->bitmask[includegroup];
}
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -82,6 +88,8 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
ipage_inner->reset();
if (respamiddle) ipage_middle->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = n_inner = 0;
neighptr = ipage->vget();
@ -95,6 +103,11 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over remaining atoms, owned and ghost
@ -110,7 +123,13 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (minchange = domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -174,27 +193,32 @@ void Neighbor::respa_nsq_no_newton(NeighList *list)
void Neighbor::respa_nsq_newton(NeighList *list)
{
int i,j,n,itype,jtype,itag,jtag,n_inner,n_middle,bitmask;
int imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int molecular = atom->molecular;
if (includegroup) {
nlocal = atom->nfirst;
bitmask = group->bitmask[includegroup];
}
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -225,6 +249,8 @@ void Neighbor::respa_nsq_newton(NeighList *list)
ipage_inner->reset();
if (respamiddle) ipage_middle->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = n_inner = 0;
neighptr = ipage->vget();
@ -239,6 +265,11 @@ void Neighbor::respa_nsq_newton(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over remaining atoms, owned and ghost
@ -270,7 +301,13 @@ void Neighbor::respa_nsq_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (minchange = domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -334,7 +371,8 @@ void Neighbor::respa_nsq_newton(NeighList *list)
void Neighbor::respa_bin_no_newton(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,n_inner,n_middle;
int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
@ -342,20 +380,23 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -388,6 +429,8 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
ipage_inner->reset();
if (respamiddle) ipage_middle->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = n_inner = 0;
neighptr = ipage->vget();
@ -402,6 +445,11 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
ytmp = x[i][1];
ztmp = x[i][2];
ibin = coord2bin(x[i]);
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in surrounding bins in stencil including self
// only store pair if i < j
@ -422,7 +470,13 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (minchange = domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -487,7 +541,8 @@ void Neighbor::respa_bin_no_newton(NeighList *list)
void Neighbor::respa_bin_newton(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,n_inner,n_middle;
int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
@ -495,20 +550,23 @@ void Neighbor::respa_bin_newton(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -541,6 +599,8 @@ void Neighbor::respa_bin_newton(NeighList *list)
ipage_inner->reset();
if (respamiddle) ipage_middle->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = n_inner = 0;
neighptr = ipage->vget();
@ -554,6 +614,11 @@ void Neighbor::respa_bin_newton(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// 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
@ -578,7 +643,13 @@ void Neighbor::respa_bin_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (minchange = domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -616,7 +687,13 @@ void Neighbor::respa_bin_newton(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (minchange = domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
@ -681,7 +758,8 @@ void Neighbor::respa_bin_newton(NeighList *list)
void Neighbor::respa_bin_newton_tri(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,n_inner,n_middle;
int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
@ -689,20 +767,23 @@ void Neighbor::respa_bin_newton_tri(NeighList *list)
bin_atoms();
// loop over each atom, storing neighbors
tagint **special = atom->special;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
int *tag = atom->tag;
int *molecule = atom->molecule;
int **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int molecular = atom->molecular;
if (includegroup) nlocal = atom->nfirst;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
int molecular = atom->molecular;
if (molecular == 2) moltemplate = 1;
else moltemplate = 0;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
@ -735,6 +816,8 @@ void Neighbor::respa_bin_newton_tri(NeighList *list)
ipage_inner->reset();
if (respamiddle) ipage_middle->reset();
// loop over each atom, storing neighbors
for (i = 0; i < nlocal; i++) {
n = n_inner = 0;
neighptr = ipage->vget();
@ -748,6 +831,11 @@ void Neighbor::respa_bin_newton_tri(NeighList *list)
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
if (moltemplate) {
imol = molindex[i];
iatom = molatom[i];
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
@ -777,7 +865,13 @@ void Neighbor::respa_bin_newton_tri(NeighList *list)
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which == 0;
if (which == 0) neighptr[n++] = j;
else if (minchange = domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;

15
src/neighbor.cpp
View File

@ -756,6 +756,7 @@ void Neighbor::init()
}
// set flags that determine which topology neighboring routines to use
// bonds,etc can only be broken for atom->molecular = 1, not 2
// SHAKE sets bonds and angles negative
// bond_quartic sets bonds to 0
// delete_bonds sets all interactions negative
@ -767,7 +768,7 @@ void Neighbor::init()
bond_off = angle_off = 1;
if (force->bond && force->bond_match("quartic")) bond_off = 1;
if (atom->avec->bonds_allow) {
if (atom->avec->bonds_allow && atom->molecular == 1) {
for (i = 0; i < atom->nlocal; i++) {
if (bond_off) break;
for (m = 0; m < atom->num_bond[i]; m++)
@ -775,7 +776,7 @@ void Neighbor::init()
}
}
if (atom->avec->angles_allow) {
if (atom->avec->angles_allow && atom->molecular == 1) {
for (i = 0; i < atom->nlocal; i++) {
if (angle_off) break;
for (m = 0; m < atom->num_angle[i]; m++)
@ -784,7 +785,7 @@ void Neighbor::init()
}
int dihedral_off = 0;
if (atom->avec->dihedrals_allow) {
if (atom->avec->dihedrals_allow && atom->molecular == 1) {
for (i = 0; i < atom->nlocal; i++) {
if (dihedral_off) break;
for (m = 0; m < atom->num_dihedral[i]; m++)
@ -793,7 +794,7 @@ void Neighbor::init()
}
int improper_off = 0;
if (atom->avec->impropers_allow) {
if (atom->avec->impropers_allow && atom->molecular == 1) {
for (i = 0; i < atom->nlocal; i++) {
if (improper_off) break;
for (m = 0; m < atom->num_improper[i]; m++)
@ -815,15 +816,19 @@ void Neighbor::init()
// set ptrs to topology build functions
if (bond_off) bond_build = &Neighbor::bond_partial;
else if (atom->molecular == 2) bond_build = &Neighbor::bond_template;
else bond_build = &Neighbor::bond_all;
if (angle_off) angle_build = &Neighbor::angle_partial;
else if (atom->molecular == 2) angle_build = &Neighbor::angle_template;
else angle_build = &Neighbor::angle_all;
if (dihedral_off) dihedral_build = &Neighbor::dihedral_partial;
else if (atom->molecular == 2) dihedral_build = &Neighbor::dihedral_template;
else dihedral_build = &Neighbor::dihedral_all;
if (improper_off) improper_build = &Neighbor::improper_partial;
else if (atom->molecular == 2) improper_build = &Neighbor::improper_template;
else improper_build = &Neighbor::improper_all;
// set topology neighbor list counts to 0
@ -1317,7 +1322,7 @@ int Neighbor::check_distance()
/* ----------------------------------------------------------------------
build all perpetual neighbor lists every few timesteps
pairwise & topology lists are created as needed
topology lists only built if topoflag = 1
topology lists only built if topoflag = 1, USER-CUDA calls with topoflag = 0
------------------------------------------------------------------------- */
void Neighbor::build(int topoflag)

4
src/neighbor.h
View File

@ -260,21 +260,25 @@ class Neighbor : protected Pointers {
BondPtr bond_build; // ptr to bond list functions
void bond_all(); // bond list with all bonds
void bond_template(); // bond list with templated bonds
void bond_partial(); // exclude certain bonds
void bond_check();
BondPtr angle_build; // ptr to angle list functions
void angle_all(); // angle list with all angles
void angle_template(); // angle list with templated bonds
void angle_partial(); // exclude certain angles
void angle_check();
BondPtr dihedral_build; // ptr to dihedral list functions
void dihedral_all(); // dihedral list with all dihedrals
void dihedral_template(); // dihedral list with templated bonds
void dihedral_partial(); // exclude certain dihedrals
void dihedral_check(int, int **);
BondPtr improper_build; // ptr to improper list functions
void improper_all(); // improper list with all impropers
void improper_template(); // improper list with templated bonds
void improper_partial(); // exclude certain impropers
// find_special: determine if atom j is in special list of atom i

171
src/read_data.cpp
View File

@ -23,6 +23,7 @@
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "molecule.h"
#include "comm.h"
#include "update.h"
#include "modify.h"
@ -220,25 +221,25 @@ void ReadData::command(int narg, char **arg)
} else if (strcmp(keyword,"Bonds") == 0) {
topoflag = bondflag = 1;
if (atom->avec->bonds_allow == 0)
if (atom->nbonds == 0)
error->all(FLERR,"Invalid data file section: Bonds");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Bonds");
bonds(firstpass);
} else if (strcmp(keyword,"Angles") == 0) {
topoflag = angleflag = 1;
if (atom->avec->angles_allow == 0)
if (atom->nangles == 0)
error->all(FLERR,"Invalid data file section: Angles");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Angles");
angles(firstpass);
} else if (strcmp(keyword,"Dihedrals") == 0) {
topoflag = dihedralflag = 1;
if (atom->avec->dihedrals_allow == 0)
if (atom->ndihedrals == 0)
error->all(FLERR,"Invalid data file section: Dihedrals");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Dihedrals");
dihedrals(firstpass);
} else if (strcmp(keyword,"Impropers") == 0) {
topoflag = improperflag = 1;
if (atom->avec->impropers_allow == 0)
if (atom->nimpropers == 0)
error->all(FLERR,"Invalid data file section: Impropers");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Impropers");
impropers(firstpass);
@ -471,10 +472,83 @@ void ReadData::command(int narg, char **arg)
// create special bond lists for molecular systems
if (atom->molecular) {
if (atom->molecular == 1) {
Special special(lmp);
special.build();
}
// for atom style template systems, count total bonds,angles,etc
if (atom->molecular == 2) {
Molecule **onemols = atom->avec->onemols;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int nlocal = atom->nlocal;
int imol,iatom;
bigint nbonds,nangles,ndihedrals,nimpropers;
nbonds = nangles = ndihedrals = nimpropers = 0;
for (int i = 0; i < nlocal; i++) {
imol = molindex[i];
iatom = molatom[i];
nbonds += onemols[imol]->num_bond[iatom];
nangles += onemols[imol]->num_angle[iatom];
ndihedrals += onemols[imol]->num_dihedral[iatom];
nimpropers += onemols[imol]->num_improper[iatom];
}
MPI_Allreduce(&nbonds,&atom->nbonds,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nangles,&atom->nangles,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&ndihedrals,&atom->ndihedrals,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nimpropers,&atom->nimpropers,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (!force->newton_bond) {
atom->nbonds /= 2;
atom->nangles /= 3;
atom->ndihedrals /= 4;
atom->nimpropers /= 4;
}
if (me == 0) {
if (atom->nbonds) {
if (screen)
fprintf(screen," " BIGINT_FORMAT " template bonds\n",atom->nbonds);
if (logfile)
fprintf(logfile," " BIGINT_FORMAT " template bonds\n",atom->nbonds);
}
if (atom->nangles) {
if (screen)
fprintf(screen," " BIGINT_FORMAT " template angles\n",
atom->nangles);
if (logfile)
fprintf(logfile," " BIGINT_FORMAT " template angles\n",
atom->nangles);
}
if (atom->ndihedrals) {
if (screen)
fprintf(screen," " BIGINT_FORMAT " template dihedrals\n",
atom->nbonds);
if (logfile)
fprintf(logfile," " BIGINT_FORMAT " template bonds\n",
atom->ndihedrals);
}
if (atom->nimpropers) {
if (screen)
fprintf(screen," " BIGINT_FORMAT " template impropers\n",
atom->nimpropers);
if (logfile)
fprintf(logfile," " BIGINT_FORMAT " template impropers\n",
atom->nimpropers);
}
}
}
// for atom style template systems
// insure nbondtypes,etc are still consistent with template molecules,
// in case data file re-defined them
if (atom->molecular == 2) atom->avec->onemols[0]->check_attributes(1);
}
/* ----------------------------------------------------------------------
@ -559,20 +633,20 @@ void ReadData::header()
// otherwise "triangles" will be matched as "angles"
} else if (strstr(line,"ellipsoids")) {
if (!avec_ellipsoid && me == 0)
error->one(FLERR,"No ellipsoids allowed with this atom style");
if (!avec_ellipsoid)
error->all(FLERR,"No ellipsoids allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nellipsoids);
} else if (strstr(line,"lines")) {
if (!avec_line && me == 0)
error->one(FLERR,"No lines allowed with this atom style");
if (!avec_line)
error->all(FLERR,"No lines allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nlines);
} else if (strstr(line,"triangles")) {
if (!avec_tri && me == 0)
error->one(FLERR,"No triangles allowed with this atom style");
if (!avec_tri)
error->all(FLERR,"No triangles allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&ntris);
} else if (strstr(line,"bodies")) {
if (!avec_body && me == 0)
error->one(FLERR,"No bodies allowed with this atom style");
if (!avec_body)
error->all(FLERR,"No bodies allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nbodies);
}
@ -620,44 +694,49 @@ void ReadData::header()
atom->nbonds < 0 || atom->nbonds >= MAXBIGINT ||
atom->nangles < 0 || atom->nangles >= MAXBIGINT ||
atom->ndihedrals < 0 || atom->ndihedrals >= MAXBIGINT ||
atom->nimpropers < 0 || atom->nimpropers >= MAXBIGINT) {
if (me == 0) error->one(FLERR,"System in data file is too big");
}
atom->nimpropers < 0 || atom->nimpropers >= MAXBIGINT)
error->all(FLERR,"System in data file is too big");
// check that exiting string is a valid section keyword
parse_keyword(1);
for (n = 0; n < NSECTIONS; n++)
if (strcmp(keyword,section_keywords[n]) == 0) break;
if (n == NSECTIONS && me == 0) {
if (n == NSECTIONS) {
char str[128];
sprintf(str,"Unknown identifier in data file: %s",keyword);
error->one(FLERR,str);
error->all(FLERR,str);
}
// error check on consistency of header values
// error checks on header values
// must be consistent with atom style and other header values
if ((atom->nbonds || atom->nbondtypes) &&
atom->avec->bonds_allow == 0 && me == 0)
error->one(FLERR,"No bonds allowed with this atom style");
atom->avec->bonds_allow == 0)
error->all(FLERR,"No bonds allowed with this atom style");
if ((atom->nangles || atom->nangletypes) &&
atom->avec->angles_allow == 0 && me == 0)
error->one(FLERR,"No angles allowed with this atom style");
atom->avec->angles_allow == 0)
error->all(FLERR,"No angles allowed with this atom style");
if ((atom->ndihedrals || atom->ndihedraltypes) &&
atom->avec->dihedrals_allow == 0 && me == 0)
error->one(FLERR,"No dihedrals allowed with this atom style");
atom->avec->dihedrals_allow == 0)
error->all(FLERR,"No dihedrals allowed with this atom style");
if ((atom->nimpropers || atom->nimpropertypes) &&
atom->avec->impropers_allow == 0 && me == 0)
error->one(FLERR,"No impropers allowed with this atom style");
atom->avec->impropers_allow == 0)
error->all(FLERR,"No impropers allowed with this atom style");
if (atom->nbonds > 0 && atom->nbondtypes <= 0 && me == 0)
error->one(FLERR,"Bonds defined but no bond types");
if (atom->nangles > 0 && atom->nangletypes <= 0 && me == 0)
error->one(FLERR,"Angles defined but no angle types");
if (atom->ndihedrals > 0 && atom->ndihedraltypes <= 0 && me == 0)
error->one(FLERR,"Dihedrals defined but no dihedral types");
if (atom->nimpropers > 0 && atom->nimpropertypes <= 0 && me == 0)
error->one(FLERR,"Impropers defined but no improper types");
if (atom->nbonds > 0 && atom->nbondtypes <= 0)
error->all(FLERR,"Bonds defined but no bond types");
if (atom->nangles > 0 && atom->nangletypes <= 0)
error->all(FLERR,"Angles defined but no angle types");
if (atom->ndihedrals > 0 && atom->ndihedraltypes <= 0)
error->all(FLERR,"Dihedrals defined but no dihedral types");
if (atom->nimpropers > 0 && atom->nimpropertypes <= 0)
error->all(FLERR,"Impropers defined but no improper types");
if (atom->molecular == 2) {
if (atom->nbonds || atom->nangles || atom->ndihedrals || atom->nimpropers)
error->all(FLERR,"No molecule topology allowed with atom style template");
}
}
/* ----------------------------------------------------------------------
@ -803,7 +882,7 @@ void ReadData::bonds(int firstpass)
if (screen) fprintf(screen," %d = max bonds/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max bonds/atom\n",maxall);
}
atom->bond_per_atom = maxall;
atom->bond_per_atom = maxall + atom->extra_bond_per_atom;
memory->destroy(count);
return;
}
@ -877,7 +956,7 @@ void ReadData::angles(int firstpass)
if (screen) fprintf(screen," %d = max angles/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max angles/atom\n",maxall);
}
atom->angle_per_atom = maxall;
atom->angle_per_atom = maxall + atom->extra_angle_per_atom;
memory->destroy(count);
return;
}
@ -951,7 +1030,7 @@ void ReadData::dihedrals(int firstpass)
if (screen) fprintf(screen," %d = max dihedrals/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max dihedrals/atom\n",maxall);
}
atom->dihedral_per_atom = maxall;
atom->dihedral_per_atom = maxall + atom->extra_dihedral_per_atom;
memory->destroy(count);
return;
}
@ -1025,7 +1104,7 @@ void ReadData::impropers(int firstpass)
if (screen) fprintf(screen," %d = max impropers/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max impropers/atom\n",maxall);
}
atom->improper_per_atom = maxall;
atom->improper_per_atom = maxall + atom->extra_improper_per_atom;
memory->destroy(count);
return;
}
@ -1095,7 +1174,7 @@ void ReadData::bonus(bigint nbonus, AtomVec *ptr, const char *type)
void ReadData::bodies(int firstpass)
{
int i,m,nchunk,nmax,ninteger,ndouble,tmp,onebody;
int i,m,nchunk,nline,nmax,ninteger,ndouble,tmp,onebody;
char *eof;
int mapflag = 0;
@ -1117,10 +1196,10 @@ void ReadData::bodies(int firstpass)
if (me == 0) {
nchunk = 0;
nlines = 0;
nline = 0;
m = 0;
while (nchunk < nmax && nlines <= CHUNK-MAXBODY) {
while (nchunk < nmax && nline <= CHUNK-MAXBODY) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(&buffer[m],"%d %d %d",&tmp,&ninteger,&ndouble);
@ -1140,7 +1219,7 @@ void ReadData::bodies(int firstpass)
}
nchunk++;
nlines += onebody+1;
nline += onebody+1;
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
@ -1334,11 +1413,11 @@ void ReadData::fix(int ifix, char *keyword)
{
int nchunk,eof;
bigint nlines = modify->fix[ifix]->read_data_skip_lines(keyword);
bigint nline = modify->fix[ifix]->read_data_skip_lines(keyword);
bigint nread = 0;
while (nread < nlines) {
nchunk = MIN(nlines-nread,CHUNK);
while (nread < nline) {
nchunk = MIN(nline-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
modify->fix[ifix]->read_data_section(keyword,nchunk,buffer);

43
src/read_restart.cpp
View File

@ -510,7 +510,7 @@ void ReadRestart::command(int narg, char **arg)
// create special bond lists for molecular systems
if (atom->molecular) {
if (atom->molecular == 1) {
Special special(lmp);
special.build();
}
@ -764,26 +764,17 @@ void ReadRestart::header(int incompatible)
domain->nonperiodic = 2;
}
// create new AtomVec class
// if style = hybrid, read additional sub-class arguments
// create new AtomVec class using any stored args
} else if (flag == ATOM_STYLE) {
char *style = read_string();
int nwords = 0;
char **words = NULL;
if (strcmp(style,"hybrid") == 0) {
nwords = read_int();
words = new char*[nwords];
for (int i = 0; i < nwords; i++) words[i] = read_string();
}
atom->create_avec(style,nwords,words);
atom->avec->read_restart_settings(fp);
for (int i = 0; i < nwords; i++) delete [] words[i];
delete [] words;
int nargcopy = read_int();
char **argcopy = new char*[nargcopy];
for (int i = 0; i < nargcopy; i++)
argcopy[i] = read_string();
atom->create_avec(style,nargcopy,argcopy);
for (int i = 0; i < nargcopy; i++) delete [] argcopy[i];
delete [] argcopy;
delete [] style;
} else if (flag == NATOMS) {
@ -933,10 +924,10 @@ void ReadRestart::file_layout()
error->all(FLERR,"Restart file is not a MPI-IO file");
if (mpiioflag) {
long *nproc_chunk_offsets;
bigint *nproc_chunk_offsets;
memory->create(nproc_chunk_offsets,nprocs,
"write_restart:nproc_chunk_offsets");
long *nproc_chunk_sizes;
bigint *nproc_chunk_sizes;
memory->create(nproc_chunk_sizes,nprocs,
"write_restart:nproc_chunk_sizes");
@ -966,7 +957,7 @@ void ReadRestart::file_layout()
}
int all_written_send_sizes_index = 0;
long current_offset = 0;
bigint current_offset = 0;
for (int i=0;i<nprocs;i++) {
nproc_chunk_offsets[i] = current_offset;
nproc_chunk_sizes[i] = 0;
@ -986,10 +977,10 @@ void ReadRestart::file_layout()
// scatter chunk sizes and offsets to all procs
MPI_Scatter(nproc_chunk_sizes, 1, MPI_LONG,
&assignedChunkSize , 1, MPI_LONG, 0,world);
MPI_Scatter(nproc_chunk_offsets, 1, MPI_LONG,
&assignedChunkOffset , 1, MPI_LONG, 0,world);
MPI_Scatter(nproc_chunk_sizes, 1, MPI_LMP_BIGINT,
&assignedChunkSize , 1, MPI_LMP_BIGINT, 0,world);
MPI_Scatter(nproc_chunk_offsets, 1, MPI_LMP_BIGINT,
&assignedChunkOffset , 1, MPI_LMP_BIGINT, 0,world);
memory->destroy(nproc_chunk_sizes);
memory->destroy(nproc_chunk_offsets);
@ -1004,7 +995,7 @@ void ReadRestart::file_layout()
if (mpiioflag) {
if (me == 0) headerOffset = ftell(fp);
MPI_Bcast(&headerOffset,1,MPI_LONG,0,world);
MPI_Bcast(&headerOffset,1,MPI_LMP_BIGINT,0,world);
}
}

2
src/read_restart.h
View File

@ -43,7 +43,7 @@ class ReadRestart : protected Pointers {
int mpiioflag; // 1 for MPIIO output, else 0
class RestartMPIIO *mpiio; // MPIIO for restart file input
int numChunksAssigned;
long assignedChunkSize;
bigint assignedChunkSize;
MPI_Offset assignedChunkOffset,headerOffset;
void file_search(char *, char *);

17
src/replicate.cpp
View File

@ -110,20 +110,11 @@ void Replicate::command(int narg, char **arg)
// old = original atom class
// atom = new replicated atom class
// if old atom style was hybrid, pass sub-style names to create_avec
Atom *old = atom;
atom = new Atom(lmp);
atom->settings(old);
int nstyles = 0;
char **keywords = NULL;
if (strcmp(old->atom_style,"hybrid") == 0) {
AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) old->avec;
nstyles = avec_hybrid->nstyles;
keywords = avec_hybrid->keywords;
}
atom->create_avec(old->atom_style,nstyles,keywords);
atom->create_avec(old->atom_style,old->avec->nargcopy,old->avec->argcopy);
// check that new system will not be too large
// new tags cannot exceed MAXTAGINT
@ -313,9 +304,10 @@ void Replicate::command(int narg, char **arg)
atom->tag[i] += atom_offset;
atom->image[i] = image;
if (atom->molecule_flag) {
if (atom->molecular) {
if (atom->molecule[i] > 0)
atom->molecule[i] += mol_offset;
if (atom->molecular == 1) {
if (atom->avec->bonds_allow)
for (j = 0; j < atom->num_bond[i]; j++)
atom->bond_atom[i][j] += atom_offset;
@ -340,6 +332,7 @@ void Replicate::command(int narg, char **arg)
atom->improper_atom4[i][j] += atom_offset;
}
}
}
} else m += static_cast<int> (buf[m]);
}
}
@ -404,7 +397,7 @@ void Replicate::command(int narg, char **arg)
// create special bond lists for molecular systems
if (atom->molecular) {
if (atom->molecular == 1) {
Special special(lmp);
special.build();
}

5
src/set.cpp
View File

@ -845,6 +845,11 @@ void Set::topology(int keyword)
{
int m,atom1,atom2,atom3,atom4;
// error check
if (atom->molecular == 2)
error->all(FLERR,"Cannot set bond topology types for atom style template");
// border swap to acquire ghost atom info
// enforce PBC before in case atoms are outside box
// init entire system since comm->exchange is done

2
src/variable.cpp
View File

@ -3541,7 +3541,7 @@ void Variable::atom_vector(char *word, Tree **tree,
error->one(FLERR,"Variable uses atom property that isn't allocated");
if (sizeof(tagint) == sizeof(smallint)) {
newtree->type = INTARRAY;
newtree->iarray = atom->molecule;
newtree->iarray = (int *) atom->molecule;
} else {
newtree->type = BIGINTARRAY;
newtree->barray = (bigint *) atom->molecule;

17
src/write_data.cpp
View File

@ -152,11 +152,11 @@ void WriteData::write(char *file)
// sum up bond,angle counts
// may be different than atom->nbonds,nangles if broken/turned-off
if (atom->nbonds || atom->nbondtypes) {
if (atom->molecular == 1 && atom->nbonds || atom->nbondtypes) {
nbonds_local = atom->avec->pack_bond(NULL);
MPI_Allreduce(&nbonds_local,&nbonds,1,MPI_LMP_BIGINT,MPI_SUM,world);
}
if (atom->nangles || atom->nangletypes) {
if (atom->molecular == 1 && atom->nangles || atom->nangletypes) {
nangles_local = atom->avec->pack_angle(NULL);
MPI_Allreduce(&nangles_local,&nangles,1,MPI_LMP_BIGINT,MPI_SUM,world);
}
@ -181,13 +181,16 @@ void WriteData::write(char *file)
}
// per atom info
// do not write molecular topology for atom_style template
if (natoms) atoms();
if (natoms) velocities();
if (atom->molecular == 1) {
if (atom->nbonds && nbonds) bonds();
if (atom->nangles && nangles) angles();
if (atom->ndihedrals) dihedrals();
if (atom->nimpropers) impropers();
}
// extra sections managed by fixes
@ -276,21 +279,19 @@ void WriteData::force_fields()
force->pair->write_data_all(fp);
}
}
if (atom->avec->bonds_allow && force->bond && force->bond->writedata) {
if (force->bond && force->bond->writedata) {
fprintf(fp,"\nBond Coeffs\n\n");
force->bond->write_data(fp);
}
if (atom->avec->angles_allow && force->angle && force->angle->writedata) {
if (force->angle && force->angle->writedata) {
fprintf(fp,"\nAngle Coeffs\n\n");
force->angle->write_data(fp);
}
if (atom->avec->dihedrals_allow && force->dihedral &&
force->dihedral->writedata) {
if (force->dihedral && force->dihedral->writedata) {
fprintf(fp,"\nDihedral Coeffs\n\n");
force->dihedral->write_data(fp);
}
if (atom->avec->impropers_allow && force->improper &&
force->improper->writedata) {
if (force->improper && force->improper->writedata) {
fprintf(fp,"\nImproper Coeffs\n\n");
force->improper->write_data(fp);
}

22
src/write_restart.cpp
View File

@ -449,25 +449,15 @@ void WriteRestart::header()
write_int(ZPERIODIC,domain->zperiodic);
write_int_vec(BOUNDARY,6,&domain->boundary[0][0]);
// atom_style must be written before atom class values
// so read_restart can create class before reading class values
// if style = hybrid, also write sub-class styles
// avec->write_restart() writes atom_style specific info
// write atom_style and its args
write_string(ATOM_STYLE,atom->atom_style);
if (strcmp(atom->atom_style,"hybrid") == 0) {
AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) atom->avec;
int nstyles = avec_hybrid->nstyles;
char **keywords = avec_hybrid->keywords;
fwrite(&nstyles,sizeof(int),1,fp);
for (int i = 0; i < nstyles; i++) {
int n = strlen(keywords[i]) + 1;
fwrite(&atom->avec->nargcopy,sizeof(int),1,fp);
for (int i = 0; i < atom->avec->nargcopy; i++) {
int n = strlen(atom->avec->argcopy[i]) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(keywords[i],sizeof(char),n,fp);
fwrite(atom->avec->argcopy[i],sizeof(char),n,fp);
}
}
atom->avec->write_restart_settings(fp);
write_bigint(NATOMS,natoms);
write_int(NTYPES,atom->ntypes);
@ -579,7 +569,7 @@ void WriteRestart::file_layout(int send_size)
if (mpiioflag) {
if (me == 0) headerOffset = ftell(fp);
MPI_Bcast(&headerOffset,1,MPI_LONG,0,world);
MPI_Bcast(&headerOffset,1,MPI_LMP_BIGINT,0,world);
}
}