From d543c0199ee3967576f2b33875480f1775e8b4a0 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 12 May 2014 22:29:09 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11973 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/delete_atoms.cpp | 133 ++++++++++++++++++++++++++++++++++++++----- src/delete_atoms.h | 1 + 2 files changed, 121 insertions(+), 13 deletions(-) diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp index 916fcadada..df2bfeae13 100644 --- a/src/delete_atoms.cpp +++ b/src/delete_atoms.cpp @@ -16,6 +16,7 @@ #include "delete_atoms.h" #include "atom.h" #include "atom_vec.h" +#include "molecule.h" #include "comm.h" #include "domain.h" #include "force.h" @@ -53,6 +54,10 @@ void DeleteAtoms::command(int narg, char **arg) // store state before delete bigint natoms_previous = atom->natoms; + bigint nbonds_previous = atom->nbonds; + bigint nangles_previous = atom->nangles; + bigint ndihedrals_previous = atom->ndihedrals; + bigint nimpropers_previous = atom->nimpropers; // delete the atoms @@ -62,6 +67,8 @@ void DeleteAtoms::command(int narg, char **arg) else if (strcmp(arg[0],"porosity") == 0) delete_porosity(narg,arg); else error->all(FLERR,"Illegal delete_atoms command"); + if (mol_flag) delete_molecule(); + // delete local atoms flagged in dlist // reset nlocal @@ -90,7 +97,7 @@ void DeleteAtoms::command(int narg, char **arg) atom->tag_extend(); } - // reset atom->natoms + // reset atom->natoms and also topology counts // reset atom->map if it exists // set nghost to 0 so old ghosts of deleted atoms won't be mapped @@ -102,17 +109,61 @@ void DeleteAtoms::command(int narg, char **arg) atom->map_set(); } - // print before and after atom count + if (mol_flag) recount_topology(); + + // print before and after atom and topology counts bigint ndelete = natoms_previous - atom->natoms; + bigint ndelete_bonds = nbonds_previous - atom->nbonds; + bigint ndelete_angles = nangles_previous - atom->nangles; + bigint ndelete_dihedrals = ndihedrals_previous - atom->ndihedrals; + bigint ndelete_impropers = nimpropers_previous - atom->nimpropers; if (comm->me == 0) { - if (screen) fprintf(screen,"Deleted " BIGINT_FORMAT - " atoms, new total = " BIGINT_FORMAT "\n", - ndelete,atom->natoms); + if (screen) { + fprintf(screen,"Deleted " BIGINT_FORMAT + " atoms, new total = " BIGINT_FORMAT "\n", + ndelete,atom->natoms); + if (mol_flag) { + if (nbonds_previous) + fprintf(screen,"Deleted " BIGINT_FORMAT + " bonds, new total = " BIGINT_FORMAT "\n", + ndelete_bonds,atom->nbonds); + if (nangles_previous) + fprintf(screen,"Deleted " BIGINT_FORMAT + " angles, new total = " BIGINT_FORMAT "\n", + ndelete_angles,atom->nangles); + if (ndihedrals_previous) + fprintf(screen,"Deleted " BIGINT_FORMAT + " dihedrals, new total = " BIGINT_FORMAT "\n", + ndelete_dihedrals,atom->ndihedrals); + if (nimpropers_previous) + fprintf(screen,"Deleted " BIGINT_FORMAT + " impropers, new total = " BIGINT_FORMAT "\n", + ndelete_impropers,atom->nimpropers); + } + } if (logfile) fprintf(logfile,"Deleted " BIGINT_FORMAT " atoms, new total = " BIGINT_FORMAT "\n", ndelete,atom->natoms); + if (mol_flag) { + if (nbonds_previous) + fprintf(logfile,"Deleted " BIGINT_FORMAT + " bonds, new total = " BIGINT_FORMAT "\n", + ndelete_bonds,atom->nbonds); + if (nangles_previous) + fprintf(logfile,"Deleted " BIGINT_FORMAT + " angles, new total = " BIGINT_FORMAT "\n", + ndelete_angles,atom->nangles); + if (ndihedrals_previous) + fprintf(logfile,"Deleted " BIGINT_FORMAT + " dihedrals, new total = " BIGINT_FORMAT "\n", + ndelete_dihedrals,atom->ndihedrals); + if (nimpropers_previous) + fprintf(logfile,"Deleted " BIGINT_FORMAT + " impropers, new total = " BIGINT_FORMAT "\n", + ndelete_impropers,atom->nimpropers); + } } } @@ -139,8 +190,6 @@ void DeleteAtoms::delete_group(int narg, char **arg) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) dlist[i] = 1; - - if (mol_flag) delete_molecule(); } /* ---------------------------------------------------------------------- @@ -167,8 +216,6 @@ void DeleteAtoms::delete_region(int narg, char **arg) for (int i = 0; i < nlocal; i++) if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) dlist[i] = 1; - - if (mol_flag) delete_molecule(); } /* ---------------------------------------------------------------------- @@ -320,8 +367,6 @@ void DeleteAtoms::delete_overlap(int narg, char **arg) break; } } - - if (mol_flag) delete_molecule(); } /* ---------------------------------------------------------------------- @@ -353,8 +398,6 @@ void DeleteAtoms::delete_porosity(int narg, char **arg) for (int i = 0; i < nlocal; i++) if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) if (random->uniform() <= porosity_fraction) dlist[i] = 1; - - if (mol_flag) delete_molecule(); } /* ---------------------------------------------------------------------- @@ -393,6 +436,67 @@ void DeleteAtoms::delete_molecule() memory->destroy(list); } +/* ---------------------------------------------------------------------- + update bond,angle,etc counts + different for atom->molecular = 1 or 2 +------------------------------------------------------------------------- */ + +void DeleteAtoms::recount_topology() +{ + bigint nbonds = 0; + bigint nangles = 0; + bigint ndihedrals = 0; + bigint nimpropers = 0; + + if (atom->molecular == 1) { + int *num_bond = atom->num_bond; + int *num_angle = atom->num_angle; + int *num_dihedral = atom->num_dihedral; + int *num_improper = atom->num_improper; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (num_bond) nbonds += num_bond[i]; + if (num_angle) nangles += num_angle[i]; + if (num_dihedral) ndihedrals += num_dihedral[i]; + if (num_improper) nimpropers += num_improper[i]; + } + + } 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; + + 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]; + } + } + + if (atom->avec->bonds_allow) + MPI_Allreduce(&nbonds,&atom->nbonds,1,MPI_LMP_BIGINT,MPI_SUM,world); + if (atom->avec->angles_allow) + MPI_Allreduce(&nangles,&atom->nangles,1,MPI_LMP_BIGINT,MPI_SUM,world); + if (atom->avec->dihedrals_allow) + MPI_Allreduce(&ndihedrals,&atom->ndihedrals,1,MPI_LMP_BIGINT,MPI_SUM,world); + if (atom->avec->impropers_allow) + 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; + } +} + /* ---------------------------------------------------------------------- callback from comm->ring() cbuf = list of N molecule IDs, put them in hash @@ -433,6 +537,9 @@ void DeleteAtoms::options(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"mol") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal delete_atoms command"); + if (atom->molecular == 0) + error->all(FLERR,"Cannot delete_atoms mol yes for " + "non-molecular systems"); if (strcmp(arg[iarg+1],"yes") == 0) mol_flag = 1; else if (strcmp(arg[iarg+1],"no") == 0) mol_flag = 0; else error->all(FLERR,"Illegal delete_atoms command"); diff --git a/src/delete_atoms.h b/src/delete_atoms.h index a583436664..cb826f09e3 100644 --- a/src/delete_atoms.h +++ b/src/delete_atoms.h @@ -40,6 +40,7 @@ class DeleteAtoms : protected Pointers { void delete_overlap(int, char **); void delete_porosity(int, char **); void delete_molecule(); + void recount_topology(); void options(int, char **); inline int sbmask(int j) {