lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity
lammps/src/special.cpp

1340 lines
37 KiB
C++
Raw Blame History

/* ----------------------------------------------------------------------
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 <cstdio>
#include "special.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "comm.h"
#include "modify.h"
#include "fix.h"
#include "accelerator_kokkos.h"
#include "atom_masks.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define RVOUS 1 // 0 for irregular, 1 for all2all
/* ---------------------------------------------------------------------- */
Special::Special(LAMMPS *lmp) : Pointers(lmp)
{
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
onetwo = onethree = onefour = NULL;
}
/* ---------------------------------------------------------------------- */
Special::~Special()
{
memory->destroy(onetwo);
memory->destroy(onethree);
memory->destroy(onefour);
}
/* ----------------------------------------------------------------------
create 1-2, 1-3, 1-4 lists of topology neighbors
store in onetwo, onethree, onefour for each atom
store 3 counters in nspecial[i]
------------------------------------------------------------------------- */
void Special::build()
{
MPI_Barrier(world);
double time1 = MPI_Wtime();
if (me == 0 && screen) {
const double * const special_lj = force->special_lj;
const double * const special_coul = force->special_coul;
fprintf(screen,"Finding 1-2 1-3 1-4 neighbors ...\n"
" special bond factors lj: %-10g %-10g %-10g\n"
" special bond factors coul: %-10g %-10g %-10g\n",
special_lj[1],special_lj[2],special_lj[3],
special_coul[1],special_coul[2],special_coul[3]);
}
// initialize nspecial counters to 0
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
nspecial[i][0] = 0;
nspecial[i][1] = 0;
nspecial[i][2] = 0;
}
// setup atomIDs and procowner vectors in rendezvous decomposition
atom_owners();
// tally nspecial[i][0] = # of 1-2 neighbors of atom i
// create onetwo[i] = list of 1-2 neighbors for atom i
if (force->newton_bond) onetwo_build_newton();
else onetwo_build_newton_off();
// print max # of 1-2 neighbors
if (me == 0) {
if (screen) fprintf(screen," %d = max # of 1-2 neighbors\n",maxall);
if (logfile) fprintf(logfile," %d = max # of 1-2 neighbors\n",maxall);
}
// done if special_bond weights for 1-3, 1-4 are set to 1.0
if (force->special_lj[2] == 1.0 && force->special_coul[2] == 1.0 &&
force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0) {
dedup();
combine();
fix_alteration();
memory->destroy(procowner);
memory->destroy(atomIDs);
timer_output(time1);
return;
}
// tally nspecial[i][1] = # of 1-3 neighbors of atom i
// create onethree[i] = list of 1-3 neighbors for atom i
onethree_build();
// print max # of 1-3 neighbors
if (me == 0) {
if (screen) fprintf(screen," %d = max # of 1-3 neighbors\n",maxall);
if (logfile) fprintf(logfile," %d = max # of 1-3 neighbors\n",maxall);
}
// done if special_bond weights for 1-4 are set to 1.0
if (force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0) {
dedup();
if (force->special_angle) angle_trim();
combine();
fix_alteration();
memory->destroy(procowner);
memory->destroy(atomIDs);
timer_output(time1);
return;
}
// tally nspecial[i][2] = # of 1-4 neighbors of atom i
// create onefour[i] = list of 1-4 neighbors for atom i
onefour_build();
// print max # of 1-4 neighbors
if (me == 0) {
if (screen) fprintf(screen," %d = max # of 1-4 neighbors\n",maxall);
if (logfile) fprintf(logfile," %d = max # of 1-4 neighbors\n",maxall);
}
// finish processing the onetwo, onethree, onefour lists
dedup();
if (force->special_angle) angle_trim();
if (force->special_dihedral) dihedral_trim();
combine();
fix_alteration();
memory->destroy(procowner);
memory->destroy(atomIDs);
timer_output(time1);
}
/* ----------------------------------------------------------------------
setup atomIDs and procowner
------------------------------------------------------------------------- */
void Special::atom_owners()
{
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
int *proclist;
memory->create(proclist,nlocal,"special:proclist");
IDRvous *idbuf = (IDRvous *)
memory->smalloc((bigint) nlocal*sizeof(IDRvous),"special:idbuf");
// setup input buf for rendezvous comm
// one datum for each owned atom: datum = owning proc, atomID
// each proc assigned every 1/Pth atom
for (int i = 0; i < nlocal; i++) {
proclist[i] = tag[i] % nprocs;
idbuf[i].me = me;
idbuf[i].atomID = tag[i];
}
// perform rendezvous operation
char *buf;
comm->rendezvous(RVOUS,nlocal,(char *) idbuf,sizeof(IDRvous),0,proclist,
rendezvous_ids,0,buf,0,(void *) this);
memory->destroy(proclist);
memory->sfree(idbuf);
}
/* ----------------------------------------------------------------------
onetwo build when newton_bond flag on
uses rendezvous comm
------------------------------------------------------------------------- */
void Special::onetwo_build_newton()
{
int i,j,m;
tagint *tag = atom->tag;
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
// nsend = # of my datums to send
int nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_bond[i]; j++) {
m = atom->map(bond_atom[i][j]);
if (m < 0 || m >= nlocal) nsend++;
}
}
int *proclist;
memory->create(proclist,nsend,"special:proclist");
PairRvous *inbuf = (PairRvous *)
memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf");
// setup input buf to rendezvous comm
// one datum for each unowned bond partner: bond partner ID, atomID
// owning proc for each datum = bond partner ID % nprocs
nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_bond[i]; j++) {
m = atom->map(bond_atom[i][j]);
if (m >= 0 && m < nlocal) continue;
proclist[nsend] = bond_atom[i][j] % nprocs;
inbuf[nsend].atomID = bond_atom[i][j];
inbuf[nsend].partnerID = tag[i];
nsend++;
}
}
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
memory->sfree(inbuf);
// set nspecial[0] and onetwo for all owned atoms
// based on owned info plus rendezvous output info
// output datums = pairs of atoms that are 1-2 neighbors
for (i = 0; i < nlocal; i++) {
nspecial[i][0] = num_bond[i];
for (j = 0; j < num_bond[i]; j++) {
m = atom->map(bond_atom[i][j]);
if (m >= 0 && m < nlocal) nspecial[m][0]++;
}
}
for (m = 0; m < nreturn; m++) {
i = atom->map(outbuf[m].atomID);
nspecial[i][0]++;
}
int max = 0;
for (i = 0; i < nlocal; i++)
max = MAX(max,nspecial[i][0]);
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
memory->create(onetwo,nlocal,maxall,"special:onetwo");
for (i = 0; i < nlocal; i++) nspecial[i][0] = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_bond[i]; j++) {
onetwo[i][nspecial[i][0]++] = bond_atom[i][j];
m = atom->map(bond_atom[i][j]);
if (m >= 0 && m < nlocal) onetwo[m][nspecial[m][0]++] = tag[i];
}
}
for (m = 0; m < nreturn; m++) {
i = atom->map(outbuf[m].atomID);
onetwo[i][nspecial[i][0]++] = outbuf[m].partnerID;
}
memory->sfree(outbuf);
}
/* ----------------------------------------------------------------------
onetwo build with newton_bond flag off
no need for rendezvous comm
------------------------------------------------------------------------- */
void Special::onetwo_build_newton_off()
{
int i,j;
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int max = 0;
for (i = 0; i < nlocal; i++)
max = MAX(max,num_bond[i]);
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
memory->create(onetwo,nlocal,maxall,"special:onetwo");
// nsend = # of my datums to send
// include nlocal datums with owner of each atom
for (i = 0; i < nlocal; i++) {
nspecial[i][0] = num_bond[i];
for (j = 0; j < num_bond[i]; j++)
onetwo[i][j] = bond_atom[i][j];
}
}
/* ----------------------------------------------------------------------
onethree build
uses rendezvous comm
------------------------------------------------------------------------- */
void Special::onethree_build()
{
int i,j,k,m,proc;
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
// nsend = # of my datums to send
int nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < nspecial[i][0]; j++) {
m = atom->map(onetwo[i][j]);
if (m < 0 || m >= nlocal) nsend += nspecial[i][0]-1;
}
}
int *proclist;
memory->create(proclist,nsend,"special:proclist");
PairRvous *inbuf = (PairRvous *)
memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf");
// setup input buf to rendezvous comm
// datums = pairs of onetwo partners where either is unknown
// these pairs are onethree neighbors
// datum = onetwo ID, onetwo ID
// owning proc for each datum = onetwo ID % nprocs
nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < nspecial[i][0]; j++) {
m = atom->map(onetwo[i][j]);
if (m >= 0 && m < nlocal) continue;
proc = onetwo[i][j] % nprocs;
for (k = 0; k < nspecial[i][0]; k++) {
if (j == k) continue;
proclist[nsend] = proc;
inbuf[nsend].atomID = onetwo[i][j];
inbuf[nsend].partnerID = onetwo[i][k];
nsend++;
}
}
}
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
memory->sfree(inbuf);
// set nspecial[1] and onethree for all owned atoms
// based on owned info plus rendezvous output info
// output datums = pairs of atoms that are 1-3 neighbors
for (i = 0; i < nlocal; i++) {
for (j = 0; j < nspecial[i][0]; j++) {
m = atom->map(onetwo[i][j]);
if (m >= 0 && m < nlocal) nspecial[m][1] += nspecial[i][0]-1;
}
}
for (m = 0; m < nreturn; m++) {
i = atom->map(outbuf[m].atomID);
nspecial[i][1]++;
}
int max = 0;
for (i = 0; i < nlocal; i++)
max = MAX(max,nspecial[i][1]);
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
memory->create(onethree,nlocal,maxall,"special:onethree");
for (i = 0; i < nlocal; i++) nspecial[i][1] = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < nspecial[i][0]; j++) {
m = atom->map(onetwo[i][j]);
if (m < 0 || m >= nlocal) continue;
for (k = 0; k < nspecial[i][0]; k++) {
if (j == k) continue;
onethree[m][nspecial[m][1]++] = onetwo[i][k];
}
}
}
for (m = 0; m < nreturn; m++) {
i = atom->map(outbuf[m].atomID);
onethree[i][nspecial[i][1]++] = outbuf[m].partnerID;
}
memory->sfree(outbuf);
}
/* ----------------------------------------------------------------------
onefour build
uses rendezvous comm
------------------------------------------------------------------------- */
void Special::onefour_build()
{
int i,j,k,m,proc;
tagint *tag = atom->tag;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
// nsend = # of my datums to send
int nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < nspecial[i][1]; j++) {
m = atom->map(onethree[i][j]);
if (m < 0 || m >= nlocal) nsend += nspecial[i][0];
}
}
int *proclist;
memory->create(proclist,nsend,"special:proclist");
PairRvous *inbuf = (PairRvous *)
memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf");
// setup input buf to rendezvous comm
// datums = pairs of onethree and onetwo partners where onethree is unknown
// these pairs are onefour neighbors
// datum = onetwo ID, onetwo ID
// owning proc for each datum = onethree ID % nprocs
nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < nspecial[i][1]; j++) {
m = atom->map(onethree[i][j]);
if (m >= 0 && m < nlocal) continue;
proc = onethree[i][j] % nprocs;
for (k = 0; k < nspecial[i][0]; k++) {
proclist[nsend] = proc;
inbuf[nsend].atomID = onethree[i][j];
inbuf[nsend].partnerID = onetwo[i][k];
nsend++;
}
}
}
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
memory->sfree(inbuf);
// set nspecial[2] and onefour for all owned atoms
// based on owned info plus rendezvous output info
// output datums = pairs of atoms that are 1-4 neighbors
for (i = 0; i < nlocal; i++) {
for (j = 0; j < nspecial[i][1]; j++) {
m = atom->map(onethree[i][j]);
if (m >= 0 && m < nlocal) nspecial[m][2] += nspecial[i][0];
}
}
for (m = 0; m < nreturn; m++) {
i = atom->map(outbuf[m].atomID);
nspecial[i][2]++;
}
int max = 0;
for (i = 0; i < nlocal; i++)
max = MAX(max,nspecial[i][2]);
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
memory->create(onefour,nlocal,maxall,"special:onefour");
for (i = 0; i < nlocal; i++) nspecial[i][2] = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < nspecial[i][1]; j++) {
m = atom->map(onethree[i][j]);
if (m < 0 || m >= nlocal) continue;
for (k = 0; k < nspecial[i][0]; k++) {
onefour[m][nspecial[m][2]++] = onetwo[i][k];
}
}
}
for (m = 0; m < nreturn; m++) {
i = atom->map(outbuf[m].atomID);
onefour[i][nspecial[i][2]++] = outbuf[m].partnerID;
}
memory->sfree(outbuf);
}
/* ----------------------------------------------------------------------
remove duplicates within each of onetwo, onethree, onefour individually
------------------------------------------------------------------------- */
void Special::dedup()
{
int i,j;
tagint m;
// clear map so it can be used as scratch space
atom->map_clear();
// use map to cull duplicates
// exclude original atom explicitly
// adjust onetwo, onethree, onefour values to reflect removed duplicates
// must unset map for each atom
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
int unique;
for (i = 0; i < nlocal; i++) {
unique = 0;
atom->map_one(tag[i],0);
for (j = 0; j < nspecial[i][0]; j++) {
m = onetwo[i][j];
if (atom->map(m) < 0) {
onetwo[i][unique++] = m;
atom->map_one(m,0);
}
}
nspecial[i][0] = unique;
atom->map_one(tag[i],-1);
for (j = 0; j < unique; j++) atom->map_one(onetwo[i][j],-1);
}
for (i = 0; i < nlocal; i++) {
unique = 0;
atom->map_one(tag[i],0);
for (j = 0; j < nspecial[i][1]; j++) {
m = onethree[i][j];
if (atom->map(m) < 0) {
onethree[i][unique++] = m;
atom->map_one(m,0);
}
}
nspecial[i][1] = unique;
atom->map_one(tag[i],-1);
for (j = 0; j < unique; j++) atom->map_one(onethree[i][j],-1);
}
for (i = 0; i < nlocal; i++) {
unique = 0;
atom->map_one(tag[i],0);
for (j = 0; j < nspecial[i][2]; j++) {
m = onefour[i][j];
if (atom->map(m) < 0) {
onefour[i][unique++] = m;
atom->map_one(m,0);
}
}
nspecial[i][2] = unique;
atom->map_one(tag[i],-1);
for (j = 0; j < unique; j++) atom->map_one(onefour[i][j],-1);
}
// re-create map
atom->map_init(0);
atom->nghost = 0;
atom->map_set();
}
/* ----------------------------------------------------------------------
concatenate onetwo, onethree, onefour into master atom->special list
remove duplicates between 3 lists, leave dup in first list it appears in
convert nspecial[0], nspecial[1], nspecial[2] into cumulative counters
------------------------------------------------------------------------- */
void Special::combine()
{
int i,j;
tagint m;
int me;
MPI_Comm_rank(world,&me);
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
// ----------------------------------------------------
// compute culled maxspecial = max # of special neighs of any atom
// ----------------------------------------------------
// clear map so it can be used as scratch space
atom->map_clear();
// unique = # of unique nspecial neighbors of one atom
// cull duplicates using map to check for them
// exclude original atom explicitly
// must unset map for each atom
int unique;
int maxspecial = 0;
for (i = 0; i < nlocal; i++) {
unique = 0;
atom->map_one(tag[i],0);
for (j = 0; j < nspecial[i][0]; j++) {
m = onetwo[i][j];
if (atom->map(m) < 0) {
unique++;
atom->map_one(m,0);
}
}
for (j = 0; j < nspecial[i][1]; j++) {
m = onethree[i][j];
if (atom->map(m) < 0) {
unique++;
atom->map_one(m,0);
}
}
for (j = 0; j < nspecial[i][2]; j++) {
m = onefour[i][j];
if (atom->map(m) < 0) {
unique++;
atom->map_one(m,0);
}
}
maxspecial = MAX(maxspecial,unique);
atom->map_one(tag[i],-1);
for (j = 0; j < nspecial[i][0]; j++) atom->map_one(onetwo[i][j],-1);
for (j = 0; j < nspecial[i][1]; j++) atom->map_one(onethree[i][j],-1);
for (j = 0; j < nspecial[i][2]; j++) atom->map_one(onefour[i][j],-1);
}
// if atom->maxspecial has been updated before, make certain
// we do not reset it to a smaller value. Since atom->maxspecial
// is initialized to 1, this ensures that it is larger than zero.
maxspecial = MAX(atom->maxspecial,maxspecial);
// compute global maxspecial
// add in extra factor from special_bonds command
// allocate correct special array with same nmax, new maxspecial
// previously allocated one must be destroyed
// must make AtomVec class update its ptr to special
MPI_Allreduce(&maxspecial,&atom->maxspecial,1,MPI_INT,MPI_MAX,world);
atom->maxspecial += force->special_extra;
// add force->special_extra only once
force->special_extra = 0;
if (me == 0) {
if (screen)
fprintf(screen," %d = max # of special neighbors\n",atom->maxspecial);
if (logfile)
fprintf(logfile," %d = max # of special neighbors\n",atom->maxspecial);
}
if (lmp->kokkos) {
AtomKokkos* atomKK = (AtomKokkos*) atom;
atomKK->modified(Host,SPECIAL_MASK);
atomKK->sync(Device,SPECIAL_MASK);
MemoryKokkos* memoryKK = (MemoryKokkos*) memory;
memoryKK->grow_kokkos(atomKK->k_special,atom->special,
atom->nmax,atom->maxspecial,"atom:special");
atomKK->modified(Device,SPECIAL_MASK);
atomKK->sync(Host,SPECIAL_MASK);
} else {
memory->destroy(atom->special);
memory->create(atom->special,atom->nmax,atom->maxspecial,"atom:special");
}
atom->avec->grow_reset();
tagint **special = atom->special;
// ----------------------------------------------------
// fill special array with 1-2, 1-3, 1-4 neighs for each atom
// ----------------------------------------------------
// again use map to cull duplicates
// exclude original atom explicitly
// adjust nspecial[i] values to reflect removed duplicates
// nspecial[i][1] and nspecial[i][2] now become cumulative counters
for (i = 0; i < nlocal; i++) {
unique = 0;
atom->map_one(tag[i],0);
for (j = 0; j < nspecial[i][0]; j++) {
m = onetwo[i][j];
if (atom->map(m) < 0) {
special[i][unique++] = m;
atom->map_one(m,0);
}
}
nspecial[i][0] = unique;
for (j = 0; j < nspecial[i][1]; j++) {
m = onethree[i][j];
if (atom->map(m) < 0) {
special[i][unique++] = m;
atom->map_one(m,0);
}
}
nspecial[i][1] = unique;
for (j = 0; j < nspecial[i][2]; j++) {
m = onefour[i][j];
if (atom->map(m) < 0) {
special[i][unique++] = m;
atom->map_one(m,0);
}
}
nspecial[i][2] = unique;
atom->map_one(tag[i],-1);
for (j = 0; j < nspecial[i][2]; j++) atom->map_one(special[i][j],-1);
}
// re-create map
atom->map_init(0);
atom->nghost = 0;
atom->map_set();
}
/* ----------------------------------------------------------------------
trim list of 1-3 neighbors by checking defined angles
delete a 1-3 neigh if they are not end atoms of a defined angle
and if they are not 1,3 or 2,4 atoms of a defined dihedral
uses rendezvous comm
------------------------------------------------------------------------- */
void Special::angle_trim()
{
int i,j,k,m;;
int *num_angle = atom->num_angle;
int *num_dihedral = atom->num_dihedral;
tagint **angle_atom1 = atom->angle_atom1;
tagint **angle_atom2 = atom->angle_atom2;
tagint **angle_atom3 = atom->angle_atom3;
tagint **dihedral_atom1 = atom->dihedral_atom1;
tagint **dihedral_atom2 = atom->dihedral_atom2;
tagint **dihedral_atom3 = atom->dihedral_atom3;
tagint **dihedral_atom4 = atom->dihedral_atom4;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
// stats on old 1-3 neighbor counts
double onethreecount = 0.0;
for (i = 0; i < nlocal; i++) onethreecount += nspecial[i][1];
double allcount;
MPI_Allreduce(&onethreecount,&allcount,1,MPI_DOUBLE,MPI_SUM,world);
if (me == 0) {
if (screen)
fprintf(screen,
" %g = # of 1-3 neighbors before angle trim\n",allcount);
if (logfile)
fprintf(logfile,
" %g = # of 1-3 neighbors before angle trim\n",allcount);
}
// if angles or dihedrals are defined
// rendezvous angle 1-3 and dihedral 1-3,2-4 pairs
if ((num_angle && atom->nangles) || (num_dihedral && atom->ndihedrals)) {
// nsend = # of my datums to send
// latter is only for angles or dihedrlas where I own atom2 (newton bond off)
int nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_angle[i]; j++) {
if (tag[i] != angle_atom2[i][j]) continue;
m = atom->map(angle_atom1[i][j]);
if (m < 0 || m >= nlocal) nsend++;
m = atom->map(angle_atom3[i][j]);
if (m < 0 || m >= nlocal) nsend++;
}
for (j = 0; j < num_dihedral[i]; j++) {
if (tag[i] != dihedral_atom2[i][j]) continue;
m = atom->map(dihedral_atom1[i][j]);
if (m < 0 || m >= nlocal) nsend++;
m = atom->map(dihedral_atom3[i][j]);
if (m < 0 || m >= nlocal) nsend++;
m = atom->map(dihedral_atom4[i][j]);
if (m < 0 || m >= nlocal) nsend++;
}
}
int *proclist;
memory->create(proclist,nsend,"special:proclist");
PairRvous *inbuf = (PairRvous *)
memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf");
// setup input buf to rendezvous comm
// datums = pairs of onetwo partners where either is unknown
// these pairs are onethree neighbors
// datum = onetwo ID, onetwo ID
// owning proc for each datum = onetwo ID % nprocs
nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_angle[i]; j++) {
if (tag[i] != angle_atom2[i][j]) continue;
m = atom->map(angle_atom1[i][j]);
if (m < 0 || m >= nlocal) {
proclist[nsend] = angle_atom1[i][j] % nprocs;
inbuf[nsend].atomID = angle_atom1[i][j];
inbuf[nsend].partnerID = angle_atom3[i][j];
nsend++;
}
m = atom->map(angle_atom3[i][j]);
if (m < 0 || m >= nlocal) {
proclist[nsend] = angle_atom3[i][j] % nprocs;
inbuf[nsend].atomID = angle_atom3[i][j];
inbuf[nsend].partnerID = angle_atom1[i][j];
nsend++;
}
}
for (j = 0; j < num_dihedral[i]; j++) {
if (tag[i] != dihedral_atom2[i][j]) continue;
m = atom->map(dihedral_atom1[i][j]);
if (m < 0 || m >= nlocal) {
proclist[nsend] = dihedral_atom1[i][j] % nprocs;
inbuf[nsend].atomID = dihedral_atom1[i][j];
inbuf[nsend].partnerID = dihedral_atom3[i][j];
nsend++;
}
m = atom->map(dihedral_atom3[i][j]);
if (m < 0 || m >= nlocal) {
proclist[nsend] = dihedral_atom3[i][j] % nprocs;
inbuf[nsend].atomID = dihedral_atom3[i][j];
inbuf[nsend].partnerID = dihedral_atom1[i][j];
nsend++;
}
m = atom->map(dihedral_atom4[i][j]);
if (m < 0 || m >= nlocal) {
proclist[nsend] = dihedral_atom4[i][j] % nprocs;
inbuf[nsend].atomID = dihedral_atom4[i][j];
inbuf[nsend].partnerID = dihedral_atom2[i][j];
nsend++;
}
}
}
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
memory->sfree(inbuf);
// flag all onethree atoms to keep
int max = 0;
for (i = 0; i < nlocal; i++)
max = MAX(max,nspecial[i][1]);
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
int **flag;
memory->create(flag,nlocal,maxall,"special:flag");
for (i = 0; i < nlocal; i++)
for (j = 0; j < nspecial[i][1]; j++)
flag[i][j] = 0;
// reset nspecial[1] and onethree for all owned atoms based on output info
// based on owned info plus rendezvous output info
// output datums = pairs of atoms that are 1-3 neighbors
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_angle[i]; j++) {
if (tag[i] != angle_atom2[i][j]) continue;
m = atom->map(angle_atom1[i][j]);
if (m >= 0 && m < nlocal) {
for (k = 0; k < nspecial[m][1]; k++)
if (onethree[m][k] == angle_atom3[i][j]) {
flag[m][k] = 1;
break;
}
}
m = atom->map(angle_atom3[i][j]);
if (m >= 0 && m < nlocal) {
for (k = 0; k < nspecial[m][1]; k++)
if (onethree[m][k] == angle_atom1[i][j]) {
flag[m][k] = 1;
break;
}
}
}
for (j = 0; j < num_dihedral[i]; j++) {
if (tag[i] != dihedral_atom2[i][j]) continue;
m = atom->map(dihedral_atom1[i][j]);
if (m >= 0 && m < nlocal) {
for (k = 0; k < nspecial[m][1]; k++)
if (onethree[m][k] == dihedral_atom3[i][j]) {
flag[m][k] = 1;
break;
}
}
m = atom->map(dihedral_atom3[i][j]);
if (m >= 0 && m < nlocal) {
for (k = 0; k < nspecial[m][1]; k++)
if (onethree[m][k] == dihedral_atom1[i][j]) {
flag[m][k] = 1;
break;
}
}
m = atom->map(dihedral_atom4[i][j]);
if (m >= 0 && m < nlocal) {
for (k = 0; k < nspecial[m][1]; k++)
if (onethree[m][k] == dihedral_atom2[i][j]) {
flag[m][k] = 1;
break;
}
}
}
}
for (m = 0; m < nreturn; m++) {
i = atom->map(outbuf[m].atomID);
for (k = 0; k < nspecial[i][1]; k++)
if (onethree[i][k] == outbuf[m].partnerID) {
flag[i][k] = 1;
break;
}
}
memory->destroy(outbuf);
// use flag values to compress onefour list for each atom
for (i = 0; i < nlocal; i++) {
j = 0;
while (j < nspecial[i][1]) {
if (flag[i][j] == 0) {
onethree[i][j] = onethree[i][nspecial[i][1]-1];
flag[i][j] = flag[i][nspecial[i][1]-1];
nspecial[i][1]--;
} else j++;
}
}
memory->destroy(flag);
// if no angles or dihedrals are defined, delete all 1-3 neighs
} else {
for (i = 0; i < nlocal; i++) nspecial[i][1] = 0;
}
// stats on new 1-3 neighbor counts
onethreecount = 0.0;
for (i = 0; i < nlocal; i++) onethreecount += nspecial[i][1];
MPI_Allreduce(&onethreecount,&allcount,1,MPI_DOUBLE,MPI_SUM,world);
if (me == 0) {
if (screen)
fprintf(screen,
" %g = # of 1-3 neighbors after angle trim\n",allcount);
if (logfile)
fprintf(logfile,
" %g = # of 1-3 neighbors after angle trim\n",allcount);
}
}
/* ----------------------------------------------------------------------
trim list of 1-4 neighbors by checking all defined dihedrals
delete a 1-4 neigh if they are not end atoms of a defined dihedral
uses rendezvous comm
------------------------------------------------------------------------- */
void Special::dihedral_trim()
{
int i,j,k,m;
int *num_dihedral = atom->num_dihedral;
tagint **dihedral_atom1 = atom->dihedral_atom1;
tagint **dihedral_atom2 = atom->dihedral_atom2;
tagint **dihedral_atom4 = atom->dihedral_atom4;
int **nspecial = atom->nspecial;
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
// stats on old 1-4 neighbor counts
double onefourcount = 0.0;
for (i = 0; i < nlocal; i++) onefourcount += nspecial[i][2];
double allcount;
MPI_Allreduce(&onefourcount,&allcount,1,MPI_DOUBLE,MPI_SUM,world);
if (me == 0) {
if (screen)
fprintf(screen,
" %g = # of 1-4 neighbors before dihedral trim\n",allcount);
if (logfile)
fprintf(logfile,
" %g = # of 1-4 neighbors before dihedral trim\n",allcount);
}
// if dihedrals are defined, rendezvous the dihedral 1-4 pairs
if (num_dihedral && atom->ndihedrals) {
// nsend = # of my datums to send
int nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_dihedral[i]; j++) {
if (tag[i] != dihedral_atom2[i][j]) continue;
m = atom->map(dihedral_atom1[i][j]);
if (m < 0 || m >= nlocal) nsend++;
m = atom->map(dihedral_atom4[i][j]);
if (m < 0 || m >= nlocal) nsend++;
}
}
int *proclist;
memory->create(proclist,nsend,"special:proclist");
PairRvous *inbuf = (PairRvous *)
memory->smalloc((bigint) nsend*sizeof(PairRvous),"special:inbuf");
// setup input buf to rendezvous comm
// datums = pairs of onefour atom IDs in a dihedral defined for my atoms
// only dihedrals where I own atom2 (in case newton_bond off)
// datum = atom1 ID and atom4 ID
// send the datum twice, to owner of atom1 ID and atom4 ID
// owning procs for each datum = atom1 or atom4 ID % nprocs
nsend = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_dihedral[i]; j++) {
if (tag[i] != dihedral_atom2[i][j]) continue;
m = atom->map(dihedral_atom1[i][j]);
if (m < 0 || m >= nlocal) {
proclist[nsend] = dihedral_atom1[i][j] % nprocs;
inbuf[nsend].atomID = dihedral_atom1[i][j];
inbuf[nsend].partnerID = dihedral_atom4[i][j];
nsend++;
}
m = atom->map(dihedral_atom4[i][j]);
if (m < 0 || m >= nlocal) {
proclist[nsend] = dihedral_atom4[i][j] % nprocs;
inbuf[nsend].atomID = dihedral_atom4[i][j];
inbuf[nsend].partnerID = dihedral_atom1[i][j];
nsend++;
}
}
}
// perform rendezvous operation
char *buf;
int nreturn = comm->rendezvous(RVOUS,nsend,(char *) inbuf,sizeof(PairRvous),
0,proclist,
rendezvous_pairs,0,buf,sizeof(PairRvous),
(void *) this);
PairRvous *outbuf = (PairRvous *) buf;
memory->destroy(proclist);
memory->sfree(inbuf);
// flag all of my onefour IDs to keep
int max = 0;
for (i = 0; i < nlocal; i++)
max = MAX(max,nspecial[i][2]);
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
int **flag;
memory->create(flag,nlocal,maxall,"special:flag");
for (i = 0; i < nlocal; i++)
for (j = 0; j < nspecial[i][2]; j++)
flag[i][j] = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_dihedral[i]; j++) {
if (tag[i] != dihedral_atom2[i][j]) continue;
m = atom->map(dihedral_atom1[i][j]);
if (m >= 0 && m < nlocal) {
for (k = 0; k < nspecial[m][2]; k++)
if (onefour[m][k] == dihedral_atom4[i][j]) {
flag[m][k] = 1;
break;
}
}
m = atom->map(dihedral_atom4[i][j]);
if (m >= 0 && m < nlocal) {
for (k = 0; k < nspecial[m][2]; k++)
if (onefour[m][k] == dihedral_atom1[i][j]) {
flag[m][k] = 1;
break;
}
}
}
}
for (m = 0; m < nreturn; m++) {
i = atom->map(outbuf[m].atomID);
for (k = 0; k < nspecial[i][2]; k++)
if (onefour[i][k] == outbuf[m].partnerID) {
flag[i][k] = 1;
break;
}
}
memory->destroy(outbuf);
// use flag values to compress onefour list for each atom
for (i = 0; i < nlocal; i++) {
j = 0;
while (j < nspecial[i][2]) {
if (flag[i][j] == 0) {
onefour[i][j] = onefour[i][nspecial[i][2]-1];
flag[i][j] = flag[i][nspecial[i][2]-1];
nspecial[i][2]--;
} else j++;
}
}
memory->destroy(flag);
// if no dihedrals are defined, delete all 1-4 neighs
} else {
for (i = 0; i < nlocal; i++) nspecial[i][2] = 0;
}
// stats on new 1-4 neighbor counts
onefourcount = 0.0;
for (i = 0; i < nlocal; i++) onefourcount += nspecial[i][2];
MPI_Allreduce(&onefourcount,&allcount,1,MPI_DOUBLE,MPI_SUM,world);
if (me == 0) {
if (screen)
fprintf(screen,
" %g = # of 1-4 neighbors after dihedral trim\n",allcount);
if (logfile)
fprintf(logfile,
" %g = # of 1-4 neighbors after dihedral trim\n",allcount);
}
}
/* ----------------------------------------------------------------------
process data for atoms assigned to me in rendezvous decomposition
inbuf = list of N IDRvous datums
no outbuf
------------------------------------------------------------------------- */
int Special::rendezvous_ids(int n, char *inbuf,
int &flag, int *&proclist, char *&outbuf,
void *ptr)
{
Special *sptr = (Special *) ptr;
Memory *memory = sptr->memory;
int *procowner;
tagint *atomIDs;
memory->create(procowner,n,"special:procowner");
memory->create(atomIDs,n,"special:atomIDs");
IDRvous *in = (IDRvous *) inbuf;
for (int i = 0; i < n; i++) {
procowner[i] = in[i].me;
atomIDs[i] = in[i].atomID;
}
// store rendezvous data in Special class
sptr->nrvous = n;
sptr->procowner = procowner;
sptr->atomIDs = atomIDs;
// flag = 0: no second comm needed in rendezvous
flag = 0;
return 0;
}
/* ----------------------------------------------------------------------
process data for atoms assigned to me in rendezvous decomposition
inbuf = list of N PairRvous datums
outbuf = same list of N PairRvous datums, routed to different procs
------------------------------------------------------------------------- */
int Special::rendezvous_pairs(int n, char *inbuf,
int &flag, int *&proclist, char *&outbuf,
void *ptr)
{
Special *sptr = (Special *) ptr;
Atom *atom = sptr->atom;
Memory *memory = sptr->memory;
// clear atom map so it can be used here as a hash table
// faster than an STL map for large atom counts
atom->map_clear();
// hash atom IDs stored in rendezvous decomposition
int nrvous = sptr->nrvous;
tagint *atomIDs = sptr->atomIDs;
for (int i = 0; i < nrvous; i++)
atom->map_one(atomIDs[i],i);
// proclist = owner of atomID in caller decomposition
PairRvous *in = (PairRvous *) inbuf;
int *procowner = sptr->procowner;
memory->create(proclist,n,"special:proclist");
int m;
for (int i = 0; i < n; i++) {
m = atom->map(in[i].atomID);
proclist[i] = procowner[m];
}
outbuf = inbuf;
// re-create atom map
atom->map_init(0);
atom->nghost = 0;
atom->map_set();
// flag = 1: outbuf = inbuf
flag = 1;
return n;
}
/* ----------------------------------------------------------------------
allow fixes to alter special list
currently, only fix drude does this
so that both the Drude core and electron are same level of neighbor
------------------------------------------------------------------------- */
void Special::fix_alteration()
{
for (int ifix = 0; ifix < modify->nfix; ifix++)
if (modify->fix[ifix]->special_alter_flag)
modify->fix[ifix]->rebuild_special();
}
/* ----------------------------------------------------------------------
print timing output
------------------------------------------------------------------------- */
void Special::timer_output(double time1)
{
double time2 = MPI_Wtime();
if (comm->me == 0) {
if (screen) fprintf(screen," special bonds CPU = %g secs\n",time2-time1);
if (logfile) fprintf(logfile," special bonds CPU = %g secs\n",time2-time1);
}
}
Reference in New Issue View Git Blame Copy Permalink