ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added old atom and mol id version of bond_neigh which is a lot faster than topology version

Browse Source
This commit is contained in:
phankl
2022-08-17 12:56:28 +02:00
parent e14997c597
commit 6228ca0a2a
3 changed files with 274 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

268
src/MESONT/pair_mesocnt.cpp
View File

@ -148,7 +148,7 @@ void PairMesoCNT::compute(int eflag, int vflag)
// update bond neighbor list when necessary
if (update->ntimestep == neighbor->lastcall) bond_neigh();
if (update->ntimestep == neighbor->lastcall) bond_neigh_topo();
// iterate over all bonds
@ -804,10 +804,155 @@ double PairMesoCNT::init_one(int /* i */, int /* j */)
}
/* ----------------------------------------------------------------------
update bond neighbor lists
update bond neighbor lists based on atom and mol IDs
------------------------------------------------------------------------- */
void PairMesoCNT::bond_neigh()
void PairMesoCNT::bond_neigh_id()
{
int nlocal = atom->nlocal;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
int *numneigh = list->numneigh;
int *numneigh_max;
memory->create(numneigh_max, nbondlist, "pair:numneigh_max");
for (int i = 0; i < nbondlist; i++) {
int i1 = bondlist[i][0];
int i2 = bondlist[i][1];
int numneigh1, numneigh2;
// prevent ghost atom with undefined neighbors
if (i1 > nlocal - 1)
numneigh1 = 0;
else
numneigh1 = numneigh[i1];
if (i2 > nlocal - 1)
numneigh2 = 0;
else
numneigh2 = numneigh[i2];
numneigh_max[i] = numneigh1 + numneigh2 + 2;
}
// create temporary arrays for chain creation
memory->create(reduced_nlist, nbondlist, "pair:reduced_nlist");
memory->create_ragged(reduced_neighlist, nbondlist, numneigh_max, "pair:reduced_neighlist");
// reduce neighbors to common list and find longest common list size
for (int i = 0; i < nbondlist; i++) {
int i1 = bondlist[i][0];
int i2 = bondlist[i][1];
neigh_common(i1, i2, reduced_nlist[i], reduced_neighlist[i]);
// sort list according to atom-id
sort(reduced_neighlist[i], reduced_nlist[i]);
}
// resize chain arrays
memory->destroy(numchainlist);
memory->destroy(nchainlist);
memory->destroy(endlist);
memory->destroy(chainlist);
memory->grow(selfid, nbondlist, "pair:selfid");
memory->grow(selfpos, nbondlist, 2, "pair:selfpos");
// count neighbor chains per bond
memory->create(numchainlist, nbondlist, "pair:numchainlist");
int numchain_max = 0;
for (int i = 0; i < nbondlist; i++) {
numchainlist[i] = count_chains(reduced_neighlist[i], reduced_nlist[i]);
if (numchain_max < numchainlist[i]) numchain_max = numchainlist[i];
}
// count neighbor chain lengths per bond
memory->create_ragged(nchainlist, nbondlist, numchainlist, "pair:nchainlist");
for (int i = 0; i < nbondlist; i++)
chain_lengths(reduced_neighlist[i], reduced_nlist[i], nchainlist[i]);
// create connected neighbor chains and check for ends
// MEMORY: prevent zero-size array creation for chainlist
memory->create_ragged(endlist, nbondlist, numchainlist, "pair:endlist");
if (numchain_max > 0)
memory->create_ragged(chainlist, nbondlist, numchainlist, nchainlist, "pair:chainlist");
else
memory->create(chainlist, 1, 1, 1, "pair:chainlist");
int nchain_max = 0;
for (int i = 0; i < nbondlist; i++) {
int *reduced_neigh = reduced_neighlist[i];
int *end = endlist[i];
int *nchain = nchainlist[i];
int **chain = chainlist[i];
// set up connected chains and check for ends
chain_split(reduced_neigh, reduced_nlist[i], nchain, chain, end);
// find longest chain
for (int j = 0; j < numchainlist[i]; j++)
if (nchain_max < nchain[j]) nchain_max = nchain[j];
// find selfid and selfpos
tagint tag1 = atom->tag[bondlist[i][0]];
tagint tag2 = atom->tag[bondlist[i][1]];
bool found1 = false;
bool found2 = false;
for (int j = 0; j < numchainlist[i]; j++) {
for (int k = 0; k < nchain[j]; k++) {
tagint temp_tag = atom->tag[chain[j][k]];
if (tag1 == temp_tag) {
selfid[i] = j;
selfpos[i][0] = k;
found1 = true;
}
else if (tag2 == temp_tag) {
selfid[i] = j;
selfpos[i][1] = k;
found2 = true;
}
if (found1 && found2) break;
}
if (found1 && found2) break;
}
}
// resize potential arrays
memory->grow(w, nchain_max, "pair:w");
memory->grow(wnode, nchain_max, "pair:wnode");
memory->grow(dq_w, nchain_max, 3, "pair:dq_w");
memory->grow(q1_dq_w, nchain_max, 3, 3, "pair:q1_dq_w");
memory->grow(q2_dq_w, nchain_max, 3, 3, "pair:q2_dq_w");
// destroy temporary arrays for chain creation
memory->destroy(numneigh_max);
memory->destroy(reduced_neighlist);
memory->destroy(reduced_nlist);
}
/* ----------------------------------------------------------------------
update bond neighbor lists based on bond topology
------------------------------------------------------------------------- */
void PairMesoCNT::bond_neigh_topo()
{
int nlocal = atom->nlocal;
int nghost = atom->nghost;
@ -1120,6 +1265,123 @@ void PairMesoCNT::neigh_common(int i1, int i2, int &numred, int *redlist)
}
}
/* ----------------------------------------------------------------------
count neighbor chains of given bond
------------------------------------------------------------------------- */
int PairMesoCNT::count_chains(int *redlist, int numred)
{
if (numred == 0) return 0;
tagint *tag = atom->tag;
tagint *mol = atom->molecule;
int count = 1;
// split neighbor list and count chains
for (int j = 0; j < numred - 1; j++) {
int j1 = redlist[j];
int j2 = redlist[j + 1];
if (tag[j2] - tag[j1] != 1 || mol[j1] != mol[j2]) count++;
}
return count;
}
/* ----------------------------------------------------------------------
count lengths of neighbor chains of given bond
------------------------------------------------------------------------- */
void PairMesoCNT::chain_lengths(int *redlist, int numred, int *nchain)
{
if (numred == 0) return;
tagint *tag = atom->tag;
tagint *mol = atom->molecule;
int clen = 0;
int cid = 0;
// split neighbor list into connected chains
for (int j = 0; j < numred - 1; j++) {
int j1 = redlist[j];
int j2 = redlist[j + 1];
clen++;
if (tag[j2] - tag[j1] != 1 || mol[j1] != mol[j2]) {
nchain[cid++] = clen;
clen = 0;
}
}
clen++;
nchain[cid++] = clen;
}
/* ----------------------------------------------------------------------
split neighbors into chains and identify ends
------------------------------------------------------------------------- */
void PairMesoCNT::chain_split(int *redlist, int numred, int *nchain, int **chain, int *end)
{
if (numred == 0) return;
tagint *tag = atom->tag;
tagint *mol = atom->molecule;
tagint *type = atom->type;
int clen = 0;
int cid = 0;
// split neighbor list into connected chains
for (int j = 0; j < numred - 1; j++) {
int j1 = redlist[j];
int j2 = redlist[j + 1];
chain[cid][clen++] = j1;
if (tag[j2] - tag[j1] != 1 || mol[j1] != mol[j2]) {
cid++;
clen = 0;
}
}
chain[cid][clen++] = redlist[numred - 1];
cid++;
// check for chain ends
for (int j = 0; j < cid; j++) {
int clen = nchain[j];
int cstart = chain[j][0];
int cend = chain[j][clen - 1];
bool estart = match_end(type[cstart]);
bool eend = match_end(type[cend]);
if (estart && eend) end[j] = 3;
else if (estart) end[j] = 1;
else if (eend) end[j] = 2;
else end[j] = 0;
}
}
/* ----------------------------------------------------------------------
insertion sort list according to corresponding atom ID
------------------------------------------------------------------------- */
void PairMesoCNT::sort(int *list, int size)
{
int i, j, temp1, temp2;
tagint *tag = atom->tag;
for (i = 1; i < size; i++) {
j = i;
temp1 = list[j - 1];
temp2 = list[j];
while (j > 0 && tag[temp1] > tag[temp2]) {
list[j] = temp1;
list[j - 1] = temp2;
j--;
temp1 = list[j - 1];
temp2 = list[j];
}
}
}
/* ---------------------------------------------------------------------- */

9
src/MESONT/pair_mesocnt.h
View File

@ -63,9 +63,16 @@ class PairMesoCNT : public Pair {
bool match_end(int);
int count_chains(int *, int);
void allocate();
void bond_neigh();
void bond_neigh_id();
void bond_neigh_topo();
void neigh_common(int, int, int &, int *);
void chain_lengths(int *, int, int *);
void chain_split(int *, int, int *, int **, int *);
void sort(int *, int);
void read_file(const char *);
void read_data(PotentialFileReader &, double *, double &, double &, int);
void read_data(PotentialFileReader &, double **, double &, double &, double &, double &, int);

2
src/MESONT/pair_mesocnt_viscous.cpp
View File

@ -82,7 +82,7 @@ void PairMesoCNTViscous::compute(int eflag, int vflag)
// update bond neighbor list when necessary
if (update->ntimestep == neighbor->lastcall) bond_neigh();
if (update->ntimestep == neighbor->lastcall) bond_neigh_topo();
// iterate over all bonds