Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/src/delete_bonds.cpp
e530ba46f410lammm-devel

Last Change
Raw File

delete_bonds.cpp
View Options

/* ----------------------------------------------------------------------
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 <mpi.h>
#include <stdlib.h>
#include <string.h>
#include "delete_bonds.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "neighbor.h"
#include "comm.h"
#include "force.h"
#include "group.h"
#include "special.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{MULTI,ATOM,BOND,ANGLE,DIHEDRAL,IMPROPER,STATS};
/* ---------------------------------------------------------------------- */
DeleteBonds::DeleteBonds(LAMMPS *lmp) : Pointers(lmp) {}
/* ---------------------------------------------------------------------- */
void DeleteBonds::command(int narg, char **arg)
{
if (domain->box_exist == 0)
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 != 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
// comm::init needs neighbor::init needs pair::init needs kspace::init, etc
if (comm->me == 0 && screen)
fprintf(screen,"System init for delete_bonds ...\n");
lmp->init();
if (comm->me == 0 && screen) fprintf(screen,"Deleting bonds ...\n");
// identify group
int igroup = group->find(arg[0]);
if (igroup == -1) error->all(FLERR,"Cannot find delete_bonds group ID");
int groupbit = group->bitmask[igroup];
// set style and which = type value
int style = -1;
if (strcmp(arg[1],"multi") == 0) style = MULTI;
else if (strcmp(arg[1],"atom") == 0) style = ATOM;
else if (strcmp(arg[1],"bond") == 0) style = BOND;
else if (strcmp(arg[1],"angle") == 0) style = ANGLE;
else if (strcmp(arg[1],"dihedral") == 0) style = DIHEDRAL;
else if (strcmp(arg[1],"improper") == 0) style = IMPROPER;
else if (strcmp(arg[1],"stats") == 0) style = STATS;
else error->all(FLERR,"Illegal delete_bonds command");
// setup list of types (atom,bond,etc) to consider
// use force->bounds(FLERR,) to allow setting of range of types
// range can be 0 to ntypes inclusive
int *tlist = NULL;
int iarg = 2;
if (style != MULTI && style != STATS) {
if (narg < 3) error->all(FLERR,"Illegal delete_bonds command");
int n = -1;
if (style == ATOM) n = atom->ntypes;
if (style == BOND) n = atom->nbondtypes;
if (style == ANGLE) n = atom->nangletypes;
if (style == DIHEDRAL) n = atom->ndihedraltypes;
if (style == IMPROPER) n = atom->nimpropertypes;
tlist = new int[n+1];
for (int i = 0; i <= n; i++) tlist[i] = 0;
int nlo,nhi;
force->bounds(FLERR,arg[2],n,nlo,nhi,0);
for (int i = nlo; i <= nhi; i++) tlist[i] = 1;
iarg++;
}
// grab optional keywords
int any_flag = 0;
int undo_flag = 0;
int remove_flag = 0;
int special_flag = 0;
int induce_flag = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"any") == 0) any_flag = 1;
else if (strcmp(arg[iarg],"undo") == 0) undo_flag = 1;
else if (strcmp(arg[iarg],"remove") == 0) remove_flag = 1;
else if (strcmp(arg[iarg],"special") == 0) special_flag = 1;
else if (strcmp(arg[iarg],"induce") == 0) induce_flag = 1;
else error->all(FLERR,"Illegal delete_bonds command");
iarg++;
}
// border swap to insure type and mask is current for off-proc atoms
// enforce PBC before in case atoms are outside box
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
// set topology interactions either off or on
// criteria for an interaction to potentially be changed (set flag = 1)
// all atoms or any atom in interaction must be in group, based on any_flag
// for style = MULTI, all bond/angle/dihedral/improper, no other criteria
// for style = ATOM, same as MULTI, plus at least one atom is specified type
// for style = BOND/ANGLE/DIHEDRAL/IMPROPER, interaction is specified type
// for style = STATS only compute stats, flag is always 0
// if flag = 1
// set interaction type negative if undo_flag = 0
// set interaction type positive if undo_flag = 1
int *mask = atom->mask;
int *type = atom->type;
int nlocal = atom->nlocal;
int i,m,n,consider,flag,itype;
int atom1,atom2,atom3,atom4;
if (atom->avec->bonds_allow &&
(style == BOND || style == MULTI || style == ATOM)) {
int *num_bond = atom->num_bond;
int **bond_type = atom->bond_type;
for (i = 0; i < nlocal; i++) {
for (m = 0; m < num_bond[i]; m++) {
atom1 = atom->map(atom->bond_atom[i][m]);
if (atom1 == -1) error->one(FLERR,"Bond atom missing in delete_bonds");
consider = 0;
if (!any_flag && mask[i] & groupbit && mask[atom1] & groupbit)
consider = 1;
if (any_flag && (mask[i] & groupbit || mask[atom1] & groupbit))
consider = 1;
if (consider) {
flag = 0;
if (style == MULTI) flag = 1;
else if (style == ATOM) {
if (tlist[type[i]] || tlist[type[atom1]]) flag = 1;
} else if (style == BOND) {
itype = abs(bond_type[i][m]);
if (tlist[itype]) flag = 1;
}
if (flag) {
if (undo_flag == 0 && bond_type[i][m] > 0)
bond_type[i][m] = -bond_type[i][m];
if (undo_flag == 1 && bond_type[i][m] < 0)
bond_type[i][m] = -bond_type[i][m];
}
}
}
}
}
if (atom->avec->angles_allow &&
(style == ANGLE || style == MULTI || style == ATOM)) {
int *num_angle = atom->num_angle;
int **angle_type = atom->angle_type;
for (i = 0; i < nlocal; i++) {
for (m = 0; m < num_angle[i]; m++) {
atom1 = atom->map(atom->angle_atom1[i][m]);
atom2 = atom->map(atom->angle_atom2[i][m]);
atom3 = atom->map(atom->angle_atom3[i][m]);
if (atom1 == -1 || atom2 == -1 || atom3 == -1)
error->one(FLERR,"Angle atom missing in delete_bonds");
consider = 0;
if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
mask[atom3] & groupbit) consider = 1;
if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
mask[atom3] & groupbit)) consider = 1;
if (consider) {
flag = 0;
if (style == MULTI) flag = 1;
else if (style == ATOM) {
if (tlist[type[atom1]] || tlist[type[atom2]] ||
tlist[type[atom3]]) flag = 1;
} else if (style == ANGLE) {
itype = abs(angle_type[i][m]);
if (tlist[itype]) flag = 1;
}
if (flag) {
if (undo_flag == 0 && angle_type[i][m] > 0)
angle_type[i][m] = -angle_type[i][m];
if (undo_flag == 1 && angle_type[i][m] < 0)
angle_type[i][m] = -angle_type[i][m];
}
}
}
}
}
if (atom->avec->dihedrals_allow &&
(style == DIHEDRAL || style == MULTI || style == ATOM)) {
int *num_dihedral = atom->num_dihedral;
int **dihedral_type = atom->dihedral_type;
for (i = 0; i < nlocal; i++) {
for (m = 0; m < num_dihedral[i]; m++) {
atom1 = atom->map(atom->dihedral_atom1[i][m]);
atom2 = atom->map(atom->dihedral_atom2[i][m]);
atom3 = atom->map(atom->dihedral_atom3[i][m]);
atom4 = atom->map(atom->dihedral_atom4[i][m]);
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
error->one(FLERR,"Dihedral atom missing in delete_bonds");
consider = 0;
if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
mask[atom3] & groupbit && mask[atom4] & groupbit) consider = 1;
if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
mask[atom3] & groupbit || mask[atom4] & groupbit))
consider = 1;
if (consider) {
flag = 0;
if (style == MULTI) flag = 1;
else if (style == ATOM) {
if (tlist[type[atom1]] || tlist[type[atom2]] ||
tlist[type[atom3]] || tlist[type[atom4]]) flag = 1;
} else if (style == DIHEDRAL) {
itype = abs(dihedral_type[i][m]);
if (tlist[itype]) flag = 1;
}
if (flag) {
if (undo_flag == 0 && dihedral_type[i][m] > 0)
dihedral_type[i][m] = -dihedral_type[i][m];
if (undo_flag == 1 && dihedral_type[i][m] < 0)
dihedral_type[i][m] = -dihedral_type[i][m];
}
}
}
}
}
if (atom->avec->impropers_allow &&
(style == IMPROPER || style == MULTI || style == ATOM)) {
int *num_improper = atom->num_improper;
int **improper_type = atom->improper_type;
for (i = 0; i < nlocal; i++) {
for (m = 0; m < num_improper[i]; m++) {
atom1 = atom->map(atom->improper_atom1[i][m]);
atom2 = atom->map(atom->improper_atom2[i][m]);
atom3 = atom->map(atom->improper_atom3[i][m]);
atom4 = atom->map(atom->improper_atom4[i][m]);
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
error->one(FLERR,"Improper atom missing in delete_bonds");
consider = 0;
if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
mask[atom3] & groupbit && mask[atom4] & groupbit) consider = 1;
if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
mask[atom3] & groupbit || mask[atom4] & groupbit))
consider = 1;
if (consider) {
flag = 0;
if (style == MULTI) flag = 1;
else if (style == ATOM) {
if (tlist[type[atom1]] || tlist[type[atom2]] ||
tlist[type[atom3]] || tlist[type[atom4]]) flag = 1;
} else if (style == IMPROPER) {
itype = abs(improper_type[i][m]);
if (tlist[itype]) flag = 1;
}
if (flag) {
if (undo_flag == 0 && improper_type[i][m] > 0)
improper_type[i][m] = -improper_type[i][m];
if (undo_flag == 1 && improper_type[i][m] < 0)
improper_type[i][m] = -improper_type[i][m];
}
}
}
}
}
delete [] tlist;
// induce turn off of angles, dihedral, impropers due to turned off bonds
// induce turn off of dihedrals due to turned off angles
// all atoms or any atom in interaction must be in group, based on any_flag
if (induce_flag) {
// circulate list of turned off bonds around ring of procs
// circulate list of turned off angles around ring of procs
}
// remove interactions if requested
// all atoms or any atom in interaction must be in group, based on any_flag
if (remove_flag) {
if (atom->avec->bonds_allow) {
for (i = 0; i < nlocal; i++) {
m = 0;
while (m < atom->num_bond[i]) {
if (atom->bond_type[i][m] <= 0) {
atom1 = atom->map(atom->bond_atom[i][m]);
flag = 0;
if (!any_flag && mask[i] & groupbit && mask[atom1] & groupbit)
flag = 1;
if (any_flag && (mask[i] & groupbit || mask[atom1] & groupbit))
flag = 1;
if (flag) {
n = atom->num_bond[i];
atom->bond_type[i][m] = atom->bond_type[i][n-1];
atom->bond_atom[i][m] = atom->bond_atom[i][n-1];
atom->num_bond[i]--;
} else m++;
} else m++;
}
}
}
if (atom->avec->angles_allow) {
for (i = 0; i < nlocal; i++) {
m = 0;
while (m < atom->num_angle[i]) {
if (atom->angle_type[i][m] <= 0) {
atom1 = atom->map(atom->angle_atom1[i][m]);
atom2 = atom->map(atom->angle_atom2[i][m]);
atom3 = atom->map(atom->angle_atom3[i][m]);
flag = 0;
if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
mask[atom3] & groupbit) flag = 1;
if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
mask[atom3] & groupbit)) flag = 1;
if (flag) {
n = atom->num_angle[i];
atom->angle_type[i][m] = atom->angle_type[i][n-1];
atom->angle_atom1[i][m] = atom->angle_atom1[i][n-1];
atom->angle_atom2[i][m] = atom->angle_atom2[i][n-1];
atom->angle_atom3[i][m] = atom->angle_atom3[i][n-1];
atom->num_angle[i]--;
} else m++;
} else m++;
}
}
}
if (atom->avec->dihedrals_allow) {
for (i = 0; i < nlocal; i++) {
m = 0;
while (m < atom->num_dihedral[i]) {
if (atom->dihedral_type[i][m] <= 0) {
atom1 = atom->map(atom->dihedral_atom1[i][m]);
atom2 = atom->map(atom->dihedral_atom2[i][m]);
atom3 = atom->map(atom->dihedral_atom3[i][m]);
atom4 = atom->map(atom->dihedral_atom4[i][m]);
flag = 0;
if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
mask[atom3] & groupbit && mask[atom4] & groupbit) flag = 1;
if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
mask[atom3] & groupbit || mask[atom4] & groupbit))
flag = 1;
if (flag) {
n = atom->num_dihedral[i];
atom->dihedral_type[i][m] = atom->dihedral_type[i][n-1];
atom->dihedral_atom1[i][m] = atom->dihedral_atom1[i][n-1];
atom->dihedral_atom2[i][m] = atom->dihedral_atom2[i][n-1];
atom->dihedral_atom3[i][m] = atom->dihedral_atom3[i][n-1];
atom->dihedral_atom4[i][m] = atom->dihedral_atom4[i][n-1];
atom->num_dihedral[i]--;
} else m++;
} else m++;
}
}
}
if (atom->avec->impropers_allow) {
for (i = 0; i < nlocal; i++) {
m = 0;
while (m < atom->num_improper[i]) {
if (atom->improper_type[i][m] <= 0) {
atom1 = atom->map(atom->improper_atom1[i][m]);
atom2 = atom->map(atom->improper_atom2[i][m]);
atom3 = atom->map(atom->improper_atom3[i][m]);
atom4 = atom->map(atom->improper_atom4[i][m]);
flag = 0;
if (!any_flag && mask[atom1] & groupbit && mask[atom2] & groupbit &&
mask[atom3] & groupbit && mask[atom4] & groupbit) flag = 1;
if (any_flag && (mask[atom1] & groupbit || mask[atom2] & groupbit ||
mask[atom3] & groupbit || mask[atom4] & groupbit))
flag = 1;
if (flag) {
n = atom->num_improper[i];
atom->improper_type[i][m] = atom->improper_type[i][n-1];
atom->improper_atom1[i][m] = atom->improper_atom1[i][n-1];
atom->improper_atom2[i][m] = atom->improper_atom2[i][n-1];
atom->improper_atom3[i][m] = atom->improper_atom3[i][n-1];
atom->improper_atom4[i][m] = atom->improper_atom4[i][n-1];
atom->num_improper[i]--;
} else m++;
} else m++;
}
}
}
}
// if interactions were removed, recompute global counts
if (remove_flag) {
if (atom->avec->bonds_allow) {
bigint nbonds = 0;
for (i = 0; i < nlocal; i++) nbonds += atom->num_bond[i];
MPI_Allreduce(&nbonds,&atom->nbonds,1,MPI_LMP_BIGINT,
MPI_SUM,world);
if (force->newton_bond == 0) atom->nbonds /= 2;
}
if (atom->avec->angles_allow) {
bigint nangles = 0;
for (i = 0; i < nlocal; i++) nangles += atom->num_angle[i];
MPI_Allreduce(&nangles,&atom->nangles,1,MPI_LMP_BIGINT,
MPI_SUM,world);
if (force->newton_bond == 0) atom->nangles /= 3;
}
if (atom->avec->dihedrals_allow) {
bigint ndihedrals = 0;
for (i = 0; i < nlocal; i++) ndihedrals += atom->num_dihedral[i];
MPI_Allreduce(&ndihedrals,&atom->ndihedrals,
1,MPI_LMP_BIGINT,MPI_SUM,world);
if (force->newton_bond == 0) atom->ndihedrals /= 4;
}
if (atom->avec->impropers_allow) {
bigint nimpropers = 0;
for (i = 0; i < nlocal; i++) nimpropers += atom->num_improper[i];
MPI_Allreduce(&nimpropers,&atom->nimpropers,
1,MPI_LMP_BIGINT,MPI_SUM,world);
if (force->newton_bond == 0) atom->nimpropers /= 4;
}
}
// compute and print stats
bigint tmp;
bigint bond_on,bond_off;
bigint angle_on,angle_off;
bigint dihedral_on,dihedral_off;
bigint improper_on,improper_off;
if (atom->avec->bonds_allow) {
bond_on = bond_off = 0;
for (i = 0; i < nlocal; i++)
for (m = 0; m < atom->num_bond[i]; m++)
if (atom->bond_type[i][m] > 0) bond_on++;
else bond_off++;
MPI_Allreduce(&bond_on,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
bond_on = tmp;
MPI_Allreduce(&bond_off,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
bond_off = tmp;
if (force->newton_bond == 0) {
bond_on /= 2;
bond_off /= 2;
}
}
if (atom->avec->angles_allow) {
angle_on = angle_off = 0;
for (i = 0; i < nlocal; i++)
for (m = 0; m < atom->num_angle[i]; m++)
if (atom->angle_type[i][m] > 0) angle_on++;
else angle_off++;
MPI_Allreduce(&angle_on,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
angle_on = tmp;
MPI_Allreduce(&angle_off,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
angle_off = tmp;
if (force->newton_bond == 0) {
angle_on /= 3;
angle_off /= 3;
}
}
if (atom->avec->dihedrals_allow) {
dihedral_on = dihedral_off = 0;
for (i = 0; i < nlocal; i++)
for (m = 0; m < atom->num_dihedral[i]; m++)
if (atom->dihedral_type[i][m] > 0) dihedral_on++;
else dihedral_off++;
MPI_Allreduce(&dihedral_on,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
dihedral_on = tmp;
MPI_Allreduce(&dihedral_off,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
dihedral_off = tmp;
if (force->newton_bond == 0) {
dihedral_on /= 4;
dihedral_off /= 4;
}
}
if (atom->avec->impropers_allow) {
improper_on = improper_off = 0;
for (i = 0; i < nlocal; i++)
for (m = 0; m < atom->num_improper[i]; m++)
if (atom->improper_type[i][m] > 0) improper_on++;
else improper_off++;
MPI_Allreduce(&improper_on,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
improper_on = tmp;
MPI_Allreduce(&improper_off,&tmp,1,MPI_LMP_BIGINT,MPI_SUM,world);
improper_off = tmp;
if (force->newton_bond == 0) {
improper_on /= 4;
improper_off /= 4;
}
}
if (comm->me == 0) {
if (atom->avec->bonds_allow) {
if (screen) fprintf(screen,
" " BIGINT_FORMAT " total bonds, " BIGINT_FORMAT
" turned on, " BIGINT_FORMAT " turned off\n",
atom->nbonds,bond_on,bond_off);
if (logfile) fprintf(logfile,
" " BIGINT_FORMAT " total bonds, " BIGINT_FORMAT
" turned on, " BIGINT_FORMAT " turned off\n",
atom->nbonds,bond_on,bond_off);
}
if (atom->avec->angles_allow) {
if (screen) fprintf(screen,
" " BIGINT_FORMAT " total angles, " BIGINT_FORMAT
" turned on, " BIGINT_FORMAT " turned off\n",
atom->nangles,angle_on,angle_off);
if (logfile) fprintf(logfile,
" " BIGINT_FORMAT " total angles, " BIGINT_FORMAT
" turned on, " BIGINT_FORMAT " turned off\n",
atom->nangles,angle_on,angle_off);
}
if (atom->avec->dihedrals_allow) {
if (screen) fprintf(screen,
" " BIGINT_FORMAT " total dihedrals, "
BIGINT_FORMAT " turned on, " BIGINT_FORMAT
" turned off\n",
atom->ndihedrals,dihedral_on,dihedral_off);
if (logfile) fprintf(logfile,
" " BIGINT_FORMAT " total dihedrals, "
BIGINT_FORMAT " turned on, " BIGINT_FORMAT
" turned off\n",
atom->ndihedrals,dihedral_on,dihedral_off);
}
if (atom->avec->impropers_allow) {
if (screen) fprintf(screen,
" " BIGINT_FORMAT " total impropers, "
BIGINT_FORMAT " turned on, " BIGINT_FORMAT
" turned off\n",
atom->nimpropers,improper_on,improper_off);
if (logfile) fprintf(logfile,
" " BIGINT_FORMAT " total impropers, "
BIGINT_FORMAT " turned on, " BIGINT_FORMAT
" turned off\n",
atom->nimpropers,improper_on,improper_off);
}
}
// re-compute special list if requested
if (special_flag) {
Special special(lmp);
special.build();
}
}
c4science · Help