From b46ec2206e7aa7a0410788d7f4ea85c06beb1024 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 16 Apr 2014 23:21:07 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11802
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/MC/fix_bond_break.cpp | 183 ++++++++++++++++++--------------------
src/MC/fix_bond_break.h | 3 +-
2 files changed, 90 insertions(+), 96 deletions(-)
diff --git a/src/MC/fix_bond_break.cpp b/src/MC/fix_bond_break.cpp
index febe418c88..0783bab26e 100755
--- a/src/MC/fix_bond_break.cpp
+++ b/src/MC/fix_bond_break.cpp
@@ -104,11 +104,11 @@ FixBondBreak::FixBondBreak(LAMMPS *lmp, int narg, char **arg) :
maxbreak = 0;
broken = NULL;
- maxinfluenced = 0;
- influenced = NULL;
-
// copy = special list for one atom
- // may contain 1-2 neighs of all 1-3 neighs before dedup() shrinks it
+ // size = ms^2 + ms is sufficient
+ // b/c in rebuild_special() neighs of all 1-2s are added,
+ // then a dedup(), then neighs of all 1-3s are added, then final dedup()
+ // this means intermediate size cannot exceed ms^2 + ms
int maxspecial = atom->maxspecial;
copy = new tagint[maxspecial*maxspecial + maxspecial];
@@ -131,7 +131,6 @@ FixBondBreak::~FixBondBreak()
memory->destroy(finalpartner);
memory->destroy(distsq);
memory->destroy(broken);
- memory->destroy(influenced);
delete [] copy;
}
@@ -152,6 +151,15 @@ void FixBondBreak::init()
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ // enable angle/dihedral/improper breaking if any defined
+
+ if (atom->nangles) angleflag = 1;
+ else angleflag = 0;
+ if (atom->ndihedrals) dihedralflag = 1;
+ else dihedralflag = 0;
+ if (atom->nimpropers) improperflag = 1;
+ else improperflag = 0;
+
if (force->improper) {
if (force->improper_match("class2") || force->improper_match("ring"))
error->all(FLERR,"Cannot yet use fix bond/break with this "
@@ -406,37 +414,27 @@ void FixBondBreak::check_ghosts()
yes if both atom IDs appear in atom's special list
else no
if influenced:
- check for angles/dihedrals/impropers to break due to broken bond
+ check for angles/dihedrals/impropers to break due to specific broken bonds
rebuild the atom's special list of 1-2,1-3,1-4 neighs
------------------------------------------------------------------------- */
void FixBondBreak::update_topology()
{
- int i,j,k,n,influence,found;
+ int i,j,k,n,influence,influenced,found;
tagint id1,id2;
tagint *slist;
- int angles_allow = atom->avec->angles_allow;
- int dihedrals_allow = atom->avec->dihedrals_allow;
- int impropers_allow = atom->avec->impropers_allow;
-
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
int nlocal = atom->nlocal;
- if (nlocal > maxinfluenced) {
- maxinfluenced = atom->nmax;
- memory->destroy(influenced);
- memory->create(influenced,maxinfluenced,"bond/break:influenced");
- }
-
nangles = 0;
ndihedrals = 0;
nimpropers = 0;
for (i = 0; i < nlocal; i++) {
- influenced[i] = 0;
+ influenced = 0;
slist = special[i];
for (j = 0; j < nbreak; j++) {
@@ -453,40 +451,100 @@ void FixBondBreak::update_topology()
if (found == 2) influence = 1;
}
if (!influence) continue;
- influenced[i] = 1;
+ influenced = 1;
- if (angles_allow) break_angles(i,id1,id2);
- if (dihedrals_allow) break_dihedrals(i,id1,id2);
- if (impropers_allow) break_impropers(i,id1,id2);
+ if (angleflag) break_angles(i,id1,id2);
+ if (dihedralflag) break_dihedrals(i,id1,id2);
+ if (improperflag) break_impropers(i,id1,id2);
}
+
+ if (influenced) rebuild_special(i);
}
int newton_bond = force->newton_bond;
int all;
- if (angles_allow) {
+ if (angleflag) {
MPI_Allreduce(&nangles,&all,1,MPI_INT,MPI_SUM,world);
if (!newton_bond) all /= 3;
atom->nangles -= all;
}
- if (dihedrals_allow) {
+ if (dihedralflag) {
MPI_Allreduce(&ndihedrals,&all,1,MPI_INT,MPI_SUM,world);
if (!newton_bond) all /= 4;
atom->ndihedrals -= all;
}
- if (impropers_allow) {
+ if (improperflag) {
MPI_Allreduce(&nimpropers,&all,1,MPI_INT,MPI_SUM,world);
if (!newton_bond) all /= 4;
atom->nimpropers -= all;
}
+}
+
+/* ----------------------------------------------------------------------
+ re-build special list of atom M
+ does not affect 1-2 neighs (already include effects of new bond)
+ affects 1-3 and 1-4 neighs due to other atom's augmented 1-2 neighs
+------------------------------------------------------------------------- */
+
+void FixBondBreak::rebuild_special(int m)
+{
+ int i,j,n,n1,cn1,cn2,cn3;
+ tagint *slist;
+
+ tagint *tag = atom->tag;
+ int **nspecial = atom->nspecial;
+ tagint **special = atom->special;
+
+ // existing 1-2 neighs of atom M
+
+ slist = special[m];
+ n1 = nspecial[m][0];
+ cn1 = 0;
+ for (i = 0; i < n1; i++)
+ copy[cn1++] = slist[i];
+
+ // new 1-3 neighs of atom M, based on 1-2 neighs of 1-2 neighs
+ // exclude self
+ // remove duplicates after adding all possible 1-3 neighs
+
+ cn2 = cn1;
+ for (i = 0; i < cn1; i++) {
+ n = atom->map(copy[i]);
+ slist = special[n];
+ n1 = nspecial[n][0];
+ for (j = 0; j < n1; j++)
+ if (slist[j] != tag[m]) copy[cn2++] = slist[j];
+ }
+
+ cn2 = dedup(cn1,cn2,copy);
+
+ // new 1-4 neighs of atom M, based on 1-2 neighs of 1-3 neighs
+ // exclude self
+ // remove duplicates after adding all possible 1-4 neighs
- for (i = 0; i < nlocal; i++)
- if (influenced[i]) rebuild_special(i);
+ cn3 = cn2;
+ for (i = cn1; i < cn2; i++) {
+ n = atom->map(copy[i]);
+ slist = special[n];
+ n1 = nspecial[n][0];
+ for (j = 0; j < n1; j++)
+ if (slist[j] != tag[m]) copy[cn3++] = slist[j];
+ }
+
+ cn3 = dedup(cn2,cn3,copy);
+
+ // store new special list with atom M
+
+ nspecial[m][0] = cn1;
+ nspecial[m][1] = cn2;
+ nspecial[m][2] = cn3;
+ memcpy(special[m],copy,cn3*sizeof(int));
}
/* ----------------------------------------------------------------------
break any angles owned by atom M that include atom IDs 1 and 2
- angle is broken if ID1-ID2 is one of 2 bonds in linear angle
+ angle is broken if ID1-ID2 is one of 2 bonds in angle (I-J,J-K)
------------------------------------------------------------------------- */
void FixBondBreak::break_angles(int m, tagint id1, tagint id2)
@@ -523,8 +581,8 @@ void FixBondBreak::break_angles(int m, tagint id1, tagint id2)
}
/* ----------------------------------------------------------------------
- break any dihedrals owned by atom M that include atom IDs 1,2,3
- dihedral is broken if ID1-ID2 is one of 3 bonds in linear dihedral
+ break any dihedrals owned by atom M that include atom IDs 1 and 2
+ dihedral is broken if ID1-ID2 is one of 3 bonds in dihedral (I-J,J-K.K-L)
------------------------------------------------------------------------- */
void FixBondBreak::break_dihedrals(int m, tagint id1, tagint id2)
@@ -565,7 +623,7 @@ void FixBondBreak::break_dihedrals(int m, tagint id1, tagint id2)
}
/* ----------------------------------------------------------------------
- break any impropers owned by atom M that include atom IDs 1,2,3
+ break any impropers owned by atom M that include atom IDs 1 and 2
improper is broken if ID1-ID2 is one of 3 bonds in improper (I-J,I-K,I-L)
------------------------------------------------------------------------- */
@@ -606,68 +664,6 @@ void FixBondBreak::break_impropers(int m, tagint id1, tagint id2)
atom->num_improper[m] = num_improper;
}
-
-/* ----------------------------------------------------------------------
- re-build special list of atom M due to bond ID1-ID2 being formed
- does not affect 1-2 neighs (already include effects of new bond)
- may affect 1-3 and 1-4 bonds
-------------------------------------------------------------------------- */
-
-void FixBondBreak::rebuild_special(int m)
-{
- int i,j,n,n1,cn1,cn2,cn3;
- tagint *slist;
-
- tagint *tag = atom->tag;
- int **nspecial = atom->nspecial;
- tagint **special = atom->special;
-
- // new 1-2 neighs of atom M
-
- slist = special[m];
- n1 = nspecial[m][0];
- cn1 = 0;
- for (i = 0; i < n1; i++)
- copy[cn1++] = slist[i];
-
- // new 1-3 neighs of atom M, based on 1-2 neighs of 1-2 neighs
- // exclude self
- // remove duplicates after adding all possible 1-3 neighs
-
- cn2 = cn1;
- for (i = 0; i < cn1; i++) {
- n = atom->map(copy[i]);
- slist = special[n];
- n1 = nspecial[n][0];
- for (j = 0; j < n1; j++)
- if (slist[j] != tag[m]) copy[cn2++] = slist[j];
- }
-
- cn2 = dedup(cn1,cn2,copy);
-
- // new 1-4 neighs of atom M, based on 1-2 neighs of 1-3 neighs
- // exclude self
- // remove duplicates after adding all possible 1-4 neighs
-
- cn3 = cn2;
- for (i = cn1; i < cn2; i++) {
- n = atom->map(copy[i]);
- slist = special[n];
- n1 = nspecial[n][0];
- for (j = 0; j < n1; j++)
- if (slist[j] != tag[m]) copy[cn3++] = slist[j];
- }
-
- cn3 = dedup(cn2,cn3,copy);
-
- // store new special list with atom M
-
- nspecial[m][0] = cn1;
- nspecial[m][1] = cn2;
- nspecial[m][2] = cn3;
- memcpy(special[m],copy,cn3*sizeof(int));
-}
-
/* ----------------------------------------------------------------------
remove all ID duplicates in copy from Nstart:Nstop-1
compare to all previous values in copy
@@ -853,6 +849,5 @@ double FixBondBreak::memory_usage()
int nmax = atom->nmax;
double bytes = 2*nmax * sizeof(tagint);
bytes += nmax * sizeof(double);
- bytes += maxinfluenced * sizeof(int);
return bytes;
}
diff --git a/src/MC/fix_bond_break.h b/src/MC/fix_bond_break.h
index e60afae6ae..3b7a1c7c8b 100755
--- a/src/MC/fix_bond_break.h
+++ b/src/MC/fix_bond_break.h
@@ -44,6 +44,7 @@ class FixBondBreak : public Fix {
int me,nprocs;
int btype,seed;
double cutoff,cutsq,fraction;
+ int angleflag,dihedralflag,improperflag;
tagint lastcheck;
int breakcount,breakcounttotal;
@@ -54,8 +55,6 @@ class FixBondBreak : public Fix {
int nbreak,maxbreak;
tagint **broken;
- int maxinfluenced;
- int *influenced;
tagint *copy;
class RanMars *random;
--
GitLab