From 0ce16af78b1297b760fbe5c6edb616dbd7fa2095 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Sat, 25 Jan 2014 22:46:08 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11328 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/GPU/fix_gpu.cpp | 8 +- src/GRANULAR/fix_pour.cpp | 22 +- src/MC/fix_bond_break.cpp | 2 +- src/MC/fix_bond_create.cpp | 2 +- src/MC/fix_bond_swap.cpp | 5 + src/MC/fix_gcmc.cpp | 20 +- src/MC/fix_gcmc.h | 4 +- src/MISC/fix_deposit.cpp | 22 +- src/MISC/fix_evaporate.cpp | 82 +- src/MOLECULE/atom_vec_template.cpp | 867 +++++++++++++++++++++ src/MOLECULE/atom_vec_template.h | 88 +++ src/MOLECULE/bond_quartic.cpp | 12 +- src/MOLECULE/pair_hbond_dreiding_lj.cpp | 52 +- src/MOLECULE/pair_hbond_dreiding_morse.cpp | 52 +- src/MPIIO/dump_atom_mpiio.cpp | 2 +- src/RIGID/fix_rigid_small.cpp | 8 +- src/RIGID/fix_rigid_small.h | 2 +- src/RIGID/fix_shake.cpp | 252 +++--- src/RIGID/fix_shake.h | 10 +- src/atom.cpp | 90 ++- src/atom.h | 3 +- src/atom_vec.cpp | 29 +- src/atom_vec.h | 14 +- src/atom_vec_body.cpp | 57 +- src/atom_vec_body.h | 10 +- src/atom_vec_hybrid.cpp | 27 +- src/atom_vec_hybrid.h | 4 +- src/atom_vec_line.h | 2 +- src/atom_vec_tri.h | 2 +- src/compute_angle_local.cpp | 47 +- src/compute_bond_local.cpp | 39 +- src/compute_dihedral_local.cpp | 40 +- src/compute_improper_local.cpp | 40 +- src/compute_property_atom.cpp | 2 +- src/compute_property_local.cpp | 5 + src/create_atoms.cpp | 81 +- src/delete_atoms.cpp | 4 +- src/delete_atoms.h | 2 +- src/delete_bonds.cpp | 3 +- src/domain.cpp | 71 +- src/domain.h | 4 +- src/dump.cpp | 3 + src/dump_image.cpp | 35 +- src/finish.cpp | 24 +- src/input.cpp | 2 +- src/molecule.cpp | 154 ++-- src/molecule.h | 4 +- src/neigh_bond.cpp | 336 ++++++++ src/neigh_full.cpp | 167 +++- src/neigh_half_bin.cpp | 152 +++- src/neigh_half_multi.cpp | 121 ++- src/neigh_half_nsq.cpp | 105 ++- src/neigh_respa.cpp | 194 +++-- src/neighbor.cpp | 17 +- src/neighbor.h | 4 + src/read_data.cpp | 171 ++-- src/read_restart.cpp | 43 +- src/read_restart.h | 2 +- src/replicate.cpp | 63 +- src/set.cpp | 5 + src/variable.cpp | 2 +- src/write_data.cpp | 25 +- src/write_restart.cpp | 24 +- 63 files changed, 2943 insertions(+), 798 deletions(-) create mode 100644 src/MOLECULE/atom_vec_template.cpp create mode 100644 src/MOLECULE/atom_vec_template.h diff --git a/src/GPU/fix_gpu.cpp b/src/GPU/fix_gpu.cpp index 28f2906c8d..c10272e90b 100644 --- a/src/GPU/fix_gpu.cpp +++ b/src/GPU/fix_gpu.cpp @@ -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; diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 55411a735d..55e5191b0e 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -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]; diff --git a/src/MC/fix_bond_break.cpp b/src/MC/fix_bond_break.cpp index 50f718492c..53b10ddba5 100755 --- a/src/MC/fix_bond_break.cpp +++ b/src/MC/fix_bond_break.cpp @@ -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 diff --git a/src/MC/fix_bond_create.cpp b/src/MC/fix_bond_create.cpp index 015a527015..ee7cfa2767 100755 --- a/src/MC/fix_bond_create.cpp +++ b/src/MC/fix_bond_create.cpp @@ -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))) diff --git a/src/MC/fix_bond_swap.cpp b/src/MC/fix_bond_swap.cpp index 79b298391a..de5c9ebec2 100644 --- a/src/MC/fix_bond_swap.cpp +++ b/src/MC/fix_bond_swap.cpp @@ -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 diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 90202ad7eb..5fe5882b7a 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -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 (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 diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h index ee532820b8..6a91df94f8 100644 --- a/src/MC/fix_gcmc.h +++ b/src/MC/fix_gcmc.h @@ -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(); diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp index f1f16a424f..371feda868 100644 --- a/src/MISC/fix_deposit.cpp +++ b/src/MISC/fix_deposit.cpp @@ -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]; diff --git a/src/MISC/fix_evaporate.cpp b/src/MISC/fix_evaporate.cpp index 421c22e666..208c944514 100644 --- a/src/MISC/fix_evaporate.cpp +++ b/src/MISC/fix_evaporate.cpp @@ -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,57 +253,70 @@ 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 (atom->avec->bonds_allow) { - if (force->newton_bond) ndeltopo[0] += atom->num_bond[i]; - else { - for (j = 0; j < atom->num_bond[i]; j++) { - if (tag[i] < atom->bond_atom[i][j]) ndeltopo[0]++; + if (molecular == 1) { + if (atom->avec->bonds_allow) { + if (force->newton_bond) ndeltopo[0] += atom->num_bond[i]; + else { + for (j = 0; j < atom->num_bond[i]; j++) { + if (tag[i] < atom->bond_atom[i][j]) ndeltopo[0]++; + } } } - } - if (atom->avec->angles_allow) { - if (force->newton_bond) ndeltopo[1] += atom->num_angle[i]; - else { - for (j = 0; j < atom->num_angle[i]; j++) { - m = atom->map(atom->angle_atom2[i][j]); - if (m >= 0 && m < nlocal) ndeltopo[1]++; + if (atom->avec->angles_allow) { + if (force->newton_bond) ndeltopo[1] += atom->num_angle[i]; + else { + for (j = 0; j < atom->num_angle[i]; j++) { + m = atom->map(atom->angle_atom2[i][j]); + if (m >= 0 && m < nlocal) ndeltopo[1]++; + } } } - } - if (atom->avec->dihedrals_allow) { - if (force->newton_bond) ndeltopo[2] += atom->num_dihedral[i]; - else { - for (j = 0; j < atom->num_dihedral[i]; j++) { - m = atom->map(atom->dihedral_atom2[i][j]); - if (m >= 0 && m < nlocal) ndeltopo[2]++; + if (atom->avec->dihedrals_allow) { + if (force->newton_bond) ndeltopo[2] += atom->num_dihedral[i]; + else { + for (j = 0; j < atom->num_dihedral[i]; j++) { + m = atom->map(atom->dihedral_atom2[i][j]); + if (m >= 0 && m < nlocal) ndeltopo[2]++; + } } } - } - if (atom->avec->impropers_allow) { - if (force->newton_bond) ndeltopo[3] += atom->num_improper[i]; - else { - for (j = 0; j < atom->num_improper[i]; j++) { - m = atom->map(atom->improper_atom2[i][j]); - if (m >= 0 && m < nlocal) ndeltopo[3]++; + if (atom->avec->impropers_allow) { + if (force->newton_bond) ndeltopo[3] += atom->num_improper[i]; + else { + for (j = 0; j < atom->num_improper[i]; j++) { + m = atom->map(atom->improper_atom2[i][j]); + if (m >= 0 && m < nlocal) ndeltopo[3]++; + } } } + + } 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) { diff --git a/src/MOLECULE/atom_vec_template.cpp b/src/MOLECULE/atom_vec_template.cpp new file mode 100644 index 0000000000..40be41eb06 --- /dev/null +++ b/src/MOLECULE/atom_vec_template.cpp @@ -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 (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; +} diff --git a/src/MOLECULE/atom_vec_template.h b/src/MOLECULE/atom_vec_template.h new file mode 100644 index 0000000000..a5423066f4 --- /dev/null +++ b/src/MOLECULE/atom_vec_template.h @@ -0,0 +1,88 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#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. + +*/ diff --git a/src/MOLECULE/bond_quartic.cpp b/src/MOLECULE/bond_quartic.cpp index d2f7a47b20..a0d2e5b702 100755 --- a/src/MOLECULE/bond_quartic.cpp +++ b/src/MOLECULE/bond_quartic.cpp @@ -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 diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index a25e116e96..3714071301 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -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; - klist = special[i]; - knum = nspecial[i][0]; + 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]; diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp index 975bc94331..195f76a48c 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -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; - klist = special[i]; - knum = nspecial[i][0]; + 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]; diff --git a/src/MPIIO/dump_atom_mpiio.cpp b/src/MPIIO/dump_atom_mpiio.cpp index ca51a66f0d..6d60b9f3a3 100644 --- a/src/MPIIO/dump_atom_mpiio.cpp +++ b/src/MPIIO/dump_atom_mpiio.cpp @@ -44,7 +44,7 @@ DumpAtomMPIIO::DumpAtomMPIIO(LAMMPS *lmp, int narg, char **arg) : DumpAtomMPIIO::~DumpAtomMPIIO() { - if (multifile == 0) MPI_File_close(&mpifh); + if (multifile == 0) MPI_File_close(&mpifh); } /* ---------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 4d411db49a..19c6d1b403 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -1288,7 +1288,7 @@ void FixRigidSmall::create_bodies() // value = index into per-body data structure // n = # of entries in hash - hash = new std::map(); + hash = new std::map(); 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 *hash = frsptr->hash; + std::map *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 *hash = frsptr->hash; + std::map *hash = frsptr->hash; double **ctr = frsptr->ctr; tagint *idclose = frsptr->idclose; double *rsqclose = frsptr->rsqclose; diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index 392a928bee..3709c50f71 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -157,7 +157,7 @@ class FixRigidSmall : public Fix { // class data used by ring communication callbacks - std::map *hash; + std::map *hash; double **bbox; double **ctr; tagint *idclose; diff --git a/src/RIGID/fix_shake.cpp b/src/RIGID/fix_shake.cpp index fa8ae3b496..692b4ad081 100644 --- a/src/RIGID/fix_shake.cpp +++ b/src/RIGID/fix_shake.cpp @@ -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; - for (i = 0; i < nlocal; i++) max = MAX(max,nspecial[i][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,10 +722,22 @@ void FixShake::find_clusters() // set npartner and partner_tag from special arrays // ----------------------------------------------------- - 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]; + 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;; + } } // ----------------------------------------------------- @@ -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) { - tagint *tag = atom->tag; - tagint **bond_atom = atom->bond_atom; - int nbonds = atom->num_bond[i]; + int m,nbonds; + int *btype; - 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 (molecular == 1) { + tagint *tag = atom->tag; + tagint **bond_atom = atom->bond_atom; + nbonds = atom->num_bond[i]; + + 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; + } + + } 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) return m; - return -1; + + 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) { - tagint **angle_atom1 = atom->angle_atom1; - tagint **angle_atom3 = atom->angle_atom3; - int nangles = atom->num_angle[i]; + int m,nangles; + int *atype; - 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 (molecular == 1) { + tagint **angle_atom1 = atom->angle_atom1; + tagint **angle_atom3 = atom->angle_atom3; + nangles = atom->num_angle[i]; + + 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; + } + + } 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) return m; - return -1; + + 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; } /* ---------------------------------------------------------------------- diff --git a/src/RIGID/fix_shake.h b/src/RIGID/fix_shake.h index 2077d79936..5e8847dd22 100644 --- a/src/RIGID/fix_shake.h +++ b/src/RIGID/fix_shake.h @@ -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 diff --git a/src/atom.cpp b/src/atom.cpp index f25a8af453..86ae09bdf5 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -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,59 +1340,53 @@ 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]; - 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]; - 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) { - 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]; - dihedral_atom1[ilocal][i] = onemol->dihedral_atom1[iatom][i] + offset; - dihedral_atom2[ilocal][i] = onemol->dihedral_atom2[iatom][i] + offset; - dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset; - dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset; - } - } - if (onemol->improperflag) { - 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]; - improper_atom1[ilocal][i] = onemol->improper_atom1[iatom][i] + offset; - improper_atom2[ilocal][i] = onemol->improper_atom2[iatom][i] + offset; - improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset; - improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset; - } + 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; } - // 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 (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 (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]; + dihedral_atom1[ilocal][i] = onemol->dihedral_atom1[iatom][i] + offset; + dihedral_atom2[ilocal][i] = onemol->dihedral_atom2[iatom][i] + offset; + dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset; + dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset; + } + + 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]; + improper_atom1[ilocal][i] = onemol->improper_atom1[iatom][i] + offset; + improper_atom2[ilocal][i] = onemol->improper_atom2[iatom][i] + offset; + improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset; + improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset; + } 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]; diff --git a/src/atom.h b/src/atom.h index d4c07f81fb..cc044f18d9 100644 --- a/src/atom.h +++ b/src/atom.h @@ -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; diff --git a/src/atom_vec.cpp b/src/atom_vec.cpp index fba68e808e..b381845e46 100644 --- a/src/atom_vec.cpp +++ b/src/atom_vec.cpp @@ -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"); } diff --git a/src/atom_vec.h b/src/atom_vec.h index b9683d38fb..30bfae919d 100644 --- a/src/atom_vec.h +++ b/src/atom_vec.h @@ -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; diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp index 429ac39137..53808be5b4 100644 --- a/src/atom_vec_body.cpp +++ b/src/atom_vec_body.cpp @@ -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 diff --git a/src/atom_vec_body.h b/src/atom_vec_body.h index c16be4e2c8..023208d24c 100644 --- a/src/atom_vec_body.h +++ b/src/atom_vec_body.h @@ -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,11 +97,7 @@ 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) + int intdoubleratio; // sizeof(double) / sizeof(int) MyPoolChunk *icp; MyPoolChunk *dcp; diff --git a/src/atom_vec_hybrid.cpp b/src/atom_vec_hybrid.cpp index d88f67a854..1af8ee0b98 100644 --- a/src/atom_vec_hybrid.cpp +++ b/src/atom_vec_hybrid.cpp @@ -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 diff --git a/src/atom_vec_hybrid.h b/src/atom_vec_hybrid.h index 2a5113d4d8..862ff7d6bb 100644 --- a/src/atom_vec_hybrid.h +++ b/src/atom_vec_hybrid.h @@ -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;} diff --git a/src/atom_vec_line.h b/src/atom_vec_line.h index 8aec9ea8d1..bdab1fd998 100644 --- a/src/atom_vec_line.h +++ b/src/atom_vec_line.h @@ -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; diff --git a/src/atom_vec_tri.h b/src/atom_vec_tri.h index cf4af31cc2..7fbfac3554 100644 --- a/src/atom_vec_tri.h +++ b/src/atom_vec_tri.h @@ -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; diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp index 6e2cd57322..b4f080f86c 100644 --- a/src/compute_angle_local.cpp +++ b/src/compute_angle_local.cpp @@ -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 (tag[atom2] != angle_atom2[atom2][i]) continue; - atom1 = atom->map(angle_atom1[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]); + 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; - atom3 = atom->map(angle_atom3[atom2][i]); 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; diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index ad6697e048..00aefeda3a 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -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++) { - atom2 = atom->map(bond_atom[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; } diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp index 438f665da6..a43681c8f3 100644 --- a/src/compute_dihedral_local.cpp +++ b/src/compute_dihedral_local.cpp @@ -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 (tag[atom2] != dihedral_atom2[atom2][i]) continue; - atom1 = atom->map(dihedral_atom1[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]); + atom3 = atom->map(dihedral_atom3[atom2][i]); + 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; - atom3 = atom->map(dihedral_atom3[atom2][i]); if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; - atom4 = atom->map(dihedral_atom4[atom2][i]); if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; if (flag) { diff --git a/src/compute_improper_local.cpp b/src/compute_improper_local.cpp index 4e66856df3..7726781439 100644 --- a/src/compute_improper_local.cpp +++ b/src/compute_improper_local.cpp @@ -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 (tag[atom2] != improper_atom2[atom2][i]) continue; - atom1 = atom->map(improper_atom1[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]); + atom3 = atom->map(improper_atom3[atom2][i]); + 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; - atom3 = atom->map(improper_atom3[atom2][i]); if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; - atom4 = atom->map(improper_atom4[atom2][i]); if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; if (flag) { diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp index 250451950e..c9e6b126e2 100644 --- a/src/compute_property_atom.cpp +++ b/src/compute_property_atom.cpp @@ -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; diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index b9fc3ab4d0..675ca7e58c 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -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"); diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 0f6fbe6d56..7a24c5eb09 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -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,38 +346,44 @@ 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 (onemol->bondflag) - for (int j = 0; j < num_bond[ilocal]; j++) - bond_atom[ilocal][j] += offset; - if (onemol->angleflag) - for (int j = 0; j < num_angle[ilocal]; j++) { - angle_atom1[ilocal][j] += offset; - angle_atom2[ilocal][j] += offset; - angle_atom3[ilocal][j] += offset; - } - if (onemol->dihedralflag) - for (int j = 0; j < num_dihedral[ilocal]; j++) { - dihedral_atom1[ilocal][j] += offset; - dihedral_atom2[ilocal][j] += offset; - dihedral_atom3[ilocal][j] += offset; - dihedral_atom4[ilocal][j] += offset; - } - if (onemol->improperflag) - for (int j = 0; j < num_improper[ilocal]; j++) { - improper_atom1[ilocal][j] += offset; - improper_atom2[ilocal][j] += offset; - improper_atom3[ilocal][j] += offset; - improper_atom4[ilocal][j] += offset; - } - if (onemol->specialflag) - for (int j = 0; j < nspecial[ilocal][2]; j++) - special[ilocal][j] += offset; + 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; + if (onemol->angleflag) + for (int j = 0; j < num_angle[ilocal]; j++) { + angle_atom1[ilocal][j] += offset; + angle_atom2[ilocal][j] += offset; + angle_atom3[ilocal][j] += offset; + } + if (onemol->dihedralflag) + for (int j = 0; j < num_dihedral[ilocal]; j++) { + dihedral_atom1[ilocal][j] += offset; + dihedral_atom2[ilocal][j] += offset; + dihedral_atom3[ilocal][j] += offset; + dihedral_atom4[ilocal][j] += offset; + } + if (onemol->improperflag) + for (int j = 0; j < num_improper[ilocal]; j++) { + improper_atom1[ilocal][j] += offset; + improper_atom2[ilocal][j] += offset; + improper_atom3[ilocal][j] += offset; + improper_atom4[ilocal][j] += offset; + } + 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(); } diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp index 5c353c7fe5..2db870416a 100644 --- a/src/delete_atoms.cpp +++ b/src/delete_atoms.cpp @@ -184,7 +184,7 @@ void DeleteAtoms::delete_region(int narg, char **arg) memory->create(list,n,"delete_atoms:list"); n = 0; - std::map::iterator pos; + std::map::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 *hash = cptr->hash; int nlocal = cptr->atom->nlocal; diff --git a/src/delete_atoms.h b/src/delete_atoms.h index f091bbbabc..27b8d73d65 100644 --- a/src/delete_atoms.h +++ b/src/delete_atoms.h @@ -33,7 +33,7 @@ class DeleteAtoms : protected Pointers { private: int *dlist; int compress_flag,mol_flag; - std::map *hash; + std::map *hash; void delete_group(int, char **); void delete_region(int, char **); diff --git a/src/delete_bonds.cpp b/src/delete_bonds.cpp index 9e33b472b0..e8a5809fe7 100644 --- a/src/delete_bonds.cpp +++ b/src/delete_bonds.cpp @@ -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 diff --git a/src/domain.cpp b/src/domain.cpp index c3f78ff873..6788b533a5 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -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++) { - k = atom->map(bond_atom[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++) { - if (bond_type[i][j] <= 0) continue; - k = atom->map(bond_atom[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; diff --git a/src/domain.h b/src/domain.h index 2d8ce8704b..f5c5e13f34 100644 --- a/src/domain.h +++ b/src/domain.h @@ -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; diff --git a/src/dump.cpp b/src/dump.cpp index 965aa29a75..5b5321119a 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -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; diff --git a/src/dump_image.cpp b/src/dump_image.cpp index 53f09b02a8..c7bfd868d4 100644 --- a/src/dump_image.cpp +++ b/src/dump_image.cpp @@ -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++) { - atom2 = atom->map(bond_atom[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]; } diff --git a/src/finish.cpp b/src/finish.cpp index a5872ee62c..c0363187f2 100644 --- a/src/finish.cpp +++ b/src/finish.cpp @@ -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); } diff --git a/src/input.cpp b/src/input.cpp index 5638cf089f..8da5b8c941 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -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 || diff --git a/src/molecule.cpp b/src/molecule.cpp index 27a4a1237d..b1a95247a6 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -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() @@ -1049,9 +1085,24 @@ void Molecule::allocate() if (qflag) memory->create(q,natoms,"molecule:q"); if (radiusflag) memory->create(radius,natoms,"molecule:radius"); if (rmassflag) memory->create(rmass,natoms,"molecule:rmass"); - + + // 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(num_bond,natoms,"molecule:num_bond"); 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"); diff --git a/src/molecule.h b/src/molecule.h index d5b406f54a..7518e2da15 100644 --- a/src/molecule.h +++ b/src/molecule.h @@ -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; diff --git a/src/neigh_bond.cpp b/src/neigh_bond.cpp index 17a6e55fd8..809ea19841 100644 --- a/src/neigh_bond.cpp +++ b/src/neigh_bond.cpp @@ -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; diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp index 6df4bfb512..87dfd1fa09 100644 --- a/src/neigh_full.cpp +++ b/src/neigh_full.cpp @@ -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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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; diff --git a/src/neigh_half_bin.cpp b/src/neigh_half_bin.cpp index 68f7832fa2..7b7fc3833e 100644 --- a/src/neigh_half_bin.cpp +++ b/src/neigh_half_bin.cpp @@ -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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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; diff --git a/src/neigh_half_multi.cpp b/src/neigh_half_multi.cpp index 6dc03a43ef..82eb6afbd6 100644 --- a/src/neigh_half_multi.cpp +++ b/src/neigh_half_multi.cpp @@ -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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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; diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp index efdd26a311..f90512f878 100644 --- a/src/neigh_half_nsq.cpp +++ b/src/neigh_half_nsq.cpp @@ -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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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; diff --git a/src/neigh_respa.cpp b/src/neigh_respa.cpp index b2d19b11e0..237cf84d42 100644 --- a/src/neigh_respa.cpp +++ b/src/neigh_respa.cpp @@ -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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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) { - which = find_special(special[i],nspecial[i],tag[j]); + 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; diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 15ac28378b..3a05f7e85b 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -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 @@ -912,7 +917,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq) else pb = &Neighbor::half_bin_no_newton_ghost; } else if (triclinic == 0) { pb = &Neighbor::half_bin_newton; - } else if (triclinic == 1) + } else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri; } else if (rq->newton == 1) { if (triclinic == 0) pb = &Neighbor::half_bin_newton; @@ -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) diff --git a/src/neighbor.h b/src/neighbor.h index 2bf33d22e1..8812aa5455 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -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 diff --git a/src/read_data.cpp b/src/read_data.cpp index c32976d28a..583af07dbf 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -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) + 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->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->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); diff --git a/src/read_restart.cpp b/src/read_restart.cpp index 2d6db0ec87..13795b0c43 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -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;idestroy(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); } } diff --git a/src/read_restart.h b/src/read_restart.h index d61140cd4a..da279481ef 100644 --- a/src/read_restart.h +++ b/src/read_restart.h @@ -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 *); diff --git a/src/replicate.cpp b/src/replicate.cpp index b31506d6c1..b50d676561 100644 --- a/src/replicate.cpp +++ b/src/replicate.cpp @@ -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,32 +304,34 @@ 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->avec->bonds_allow) - for (j = 0; j < atom->num_bond[i]; j++) - atom->bond_atom[i][j] += atom_offset; - if (atom->avec->angles_allow) - for (j = 0; j < atom->num_angle[i]; j++) { - atom->angle_atom1[i][j] += atom_offset; - atom->angle_atom2[i][j] += atom_offset; - atom->angle_atom3[i][j] += atom_offset; - } - if (atom->avec->dihedrals_allow) - for (j = 0; j < atom->num_dihedral[i]; j++) { - atom->dihedral_atom1[i][j] += atom_offset; - atom->dihedral_atom2[i][j] += atom_offset; - atom->dihedral_atom3[i][j] += atom_offset; - atom->dihedral_atom4[i][j] += atom_offset; - } - if (atom->avec->impropers_allow) - for (j = 0; j < atom->num_improper[i]; j++) { - atom->improper_atom1[i][j] += atom_offset; - atom->improper_atom2[i][j] += atom_offset; - atom->improper_atom3[i][j] += atom_offset; - atom->improper_atom4[i][j] += atom_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; + if (atom->avec->angles_allow) + for (j = 0; j < atom->num_angle[i]; j++) { + atom->angle_atom1[i][j] += atom_offset; + atom->angle_atom2[i][j] += atom_offset; + atom->angle_atom3[i][j] += atom_offset; + } + if (atom->avec->dihedrals_allow) + for (j = 0; j < atom->num_dihedral[i]; j++) { + atom->dihedral_atom1[i][j] += atom_offset; + atom->dihedral_atom2[i][j] += atom_offset; + atom->dihedral_atom3[i][j] += atom_offset; + atom->dihedral_atom4[i][j] += atom_offset; + } + if (atom->avec->impropers_allow) + for (j = 0; j < atom->num_improper[i]; j++) { + atom->improper_atom1[i][j] += atom_offset; + atom->improper_atom2[i][j] += atom_offset; + atom->improper_atom3[i][j] += atom_offset; + atom->improper_atom4[i][j] += atom_offset; + } + } } } else m += static_cast (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(); } diff --git a/src/set.cpp b/src/set.cpp index 108aecc239..c221fbc3e0 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -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 diff --git a/src/variable.cpp b/src/variable.cpp index c866a02ad7..065d3605d2 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -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; diff --git a/src/write_data.cpp b/src/write_data.cpp index 6c9cd49894..c3b038ee41 100644 --- a/src/write_data.cpp +++ b/src/write_data.cpp @@ -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->nbonds && nbonds) bonds(); - if (atom->nangles && nangles) angles(); - if (atom->ndihedrals) dihedrals(); - if (atom->nimpropers) impropers(); + 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); } diff --git a/src/write_restart.cpp b/src/write_restart.cpp index 2be1cfecfd..81f5eb2a6e 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -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(&n,sizeof(int),1,fp); - fwrite(keywords[i],sizeof(char),n,fp); - } + 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(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); } }