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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9473 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2013-02-13 18:20:52 +00:00
parent e6caca22d1
commit fa24d1db95
5 changed files with 441 additions and 355 deletions
Download Patch File Download Diff File Expand all files Collapse all files

214
src/fix_shake.cpp
View File

@ -38,6 +38,10 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
// allocate space for static class variable
FixShake *FixShake::fsptr;
#define BIG 1.0e20
#define MASSDELTA 0.1
@ -521,7 +525,7 @@ void FixShake::post_force(int vflag)
int m;
for (int i = 0; i < nlist; i++) {
m = list[i];
if (shake_flag[m] == 2) shake2(m);
if (shake_flag[m] == 2) shake(m);
else if (shake_flag[m] == 3) shake3(m);
else if (shake_flag[m] == 4) shake4(m);
else shake3angle(m);
@ -560,7 +564,7 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
int m;
for (int i = 0; i < nlist; i++) {
m = list[i];
if (shake_flag[m] == 2) shake2(m);
if (shake_flag[m] == 2) shake(m);
else if (shake_flag[m] == 3) shake3(m);
else if (shake_flag[m] == 4) shake4(m);
else shake3angle(m);
@ -656,7 +660,7 @@ void FixShake::find_clusters()
int max = 0;
for (i = 0; i < nlocal; i++) max = MAX(max,nspecial[i][0]);
int *npartner,*nshake;
int *npartner;
memory->create(npartner,nlocal,"shake:npartner");
memory->create(nshake,nlocal,"shake:nshake");
@ -719,10 +723,7 @@ void FixShake::find_clusters()
}
}
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
buf = new int[nbufmax];
bufcopy = new int[nbufmax];
memory->create(buf,nbuf,"shake:buf");
// fill buffer with info
@ -745,39 +746,9 @@ void FixShake::find_clusters()
}
// cycle buffer around ring of procs back to self
// when receive buffer, scan bond partner IDs for atoms I own
// if I own partner:
// fill in mask and type and massflag
// search for bond with 1st atom and fill in bondtype
messtag = 1;
for (loop = 0; loop < nprocs; loop++) {
i = 0;
while (i < size) {
m = atom->map(buf[i+1]);
if (m >= 0 && m < nlocal) {
buf[i+2] = mask[m];
buf[i+3] = type[m];
if (nmass) {
if (rmass) massone = rmass[m];
else massone = mass[type[m]];
buf[i+4] = masscheck(massone);
}
if (buf[i+5] == 0) {
n = bondfind(m,buf[i],buf[i+1]);
if (n >= 0) buf[i+5] = bond_type[m][n];
}
}
i += nper;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
fsptr = this;
comm->ring(size,sizeof(int),buf,1,ring_bonds,buf);
// store partner info returned to me
@ -793,8 +764,7 @@ void FixShake::find_clusters()
m += nper;
}
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
// error check for unfilled partner info
// if partner_type not set, is an error
@ -885,10 +855,7 @@ void FixShake::find_clusters()
}
}
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
buf = new int[nbufmax];
bufcopy = new int[nbufmax];
memory->create(buf,nbuf,"shake:buf");
// fill buffer with info
@ -905,25 +872,9 @@ void FixShake::find_clusters()
}
// cycle buffer around ring of procs back to self
// when receive buffer, scan bond partner IDs for atoms I own
// if I own partner, fill in nshake value
messtag = 2;
for (loop = 0; loop < nprocs; loop++) {
i = 0;
while (i < size) {
m = atom->map(buf[i+1]);
if (m >= 0 && m < nlocal) buf[i+2] = nshake[m];
i += 3;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
fsptr = this;
comm->ring(size,sizeof(int),buf,2,ring_nshake,buf);
// store partner info returned to me
@ -936,8 +887,7 @@ void FixShake::find_clusters()
m += 3;
}
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
// -----------------------------------------------------
// error checks
@ -1050,10 +1000,7 @@ void FixShake::find_clusters()
}
}
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
buf = new int[nbufmax];
bufcopy = new int[nbufmax];
memory->create(buf,nbuf,"shake:buf");
// fill buffer with info
@ -1079,37 +1026,11 @@ void FixShake::find_clusters()
}
// cycle buffer around ring of procs back to self
// when receive buffer, scan for ID that I own
// if I own ID, fill in shake array values
messtag = 3;
for (loop = 0; loop < nprocs; loop++) {
i = 0;
while (i < size) {
m = atom->map(buf[i]);
if (m >= 0 && m < nlocal) {
shake_flag[m] = buf[i+1];
shake_atom[m][0] = buf[i+2];
shake_atom[m][1] = buf[i+3];
shake_atom[m][2] = buf[i+4];
shake_atom[m][3] = buf[i+5];
shake_type[m][0] = buf[i+6];
shake_type[m][1] = buf[i+7];
shake_type[m][2] = buf[i+8];
}
i += 9;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
fsptr = this;
comm->ring(size,sizeof(int),buf,3,ring_shake,NULL);
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
// -----------------------------------------------------
// free local memory
@ -1196,6 +1117,99 @@ void FixShake::find_clusters()
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan bond partner IDs for atoms I own
if I own partner:
fill in mask and type and massflag
search for bond with 1st atom and fill in bondtype
------------------------------------------------------------------------- */
void FixShake::ring_bonds(int ndatum, char *cbuf)
{
Atom *atom = fsptr->atom;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *mask = atom->mask;
int **bond_type = atom->bond_type;
int *type = atom->type;
int nlocal = atom->nlocal;
int nmass = fsptr->nmass;
int *buf = (int *) cbuf;
int m,n;
double massone;
for (int i = 0; i < ndatum; i += 6) {
m = atom->map(buf[i+1]);
if (m >= 0 && m < nlocal) {
buf[i+2] = mask[m];
buf[i+3] = type[m];
if (nmass) {
if (rmass) massone = rmass[m];
else massone = mass[type[m]];
buf[i+4] = fsptr->masscheck(massone);
}
if (buf[i+5] == 0) {
n = fsptr->bondfind(m,buf[i],buf[i+1]);
if (n >= 0) buf[i+5] = bond_type[m][n];
}
}
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan bond partner IDs for atoms I own
if I own partner, fill in nshake value
------------------------------------------------------------------------- */
void FixShake::ring_nshake(int ndatum, char *cbuf)
{
Atom *atom = fsptr->atom;
int nlocal = atom->nlocal;
int *nshake = fsptr->nshake;
int *buf = (int *) cbuf;
int m;
for (int i = 0; i < ndatum; i += 3) {
m = atom->map(buf[i+1]);
if (m >= 0 && m < nlocal) buf[i+2] = nshake[m];
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan bond partner IDs for atoms I own
if I own partner, fill in nshake value
------------------------------------------------------------------------- */
void FixShake::ring_shake(int ndatum, char *cbuf)
{
Atom *atom = fsptr->atom;
int nlocal = atom->nlocal;
int *shake_flag = fsptr->shake_flag;
int **shake_atom = fsptr->shake_atom;
int **shake_type = fsptr->shake_type;
int *buf = (int *) cbuf;
int m;
for (int i = 0; i < ndatum; i += 9) {
m = atom->map(buf[i]);
if (m >= 0 && m < nlocal) {
shake_flag[m] = buf[i+1];
shake_atom[m][0] = buf[i+2];
shake_atom[m][1] = buf[i+3];
shake_atom[m][2] = buf[i+4];
shake_atom[m][3] = buf[i+5];
shake_type[m][0] = buf[i+6];
shake_type[m][1] = buf[i+7];
shake_type[m][2] = buf[i+8];
}
}
}
/* ----------------------------------------------------------------------
check if massone is within MASSDELTA of any mass in mass_list
return 1 if yes, 0 if not
@ -1252,7 +1266,7 @@ void FixShake::unconstrained_update_respa(int ilevel)
// for levels > 0 this includes more than one velocity update
// xshake = predicted position from call to this routine at level N =
// x + dt0 (v + dtN/m fN + 1/2 dt(N-1)/m f(N-1) + ... + 1/2 dt0/m f0)
// also set dtfsq = dt0*dtN so that shake2,shake3,etc can use it
// also set dtfsq = dt0*dtN so that shake,shake3,etc can use it
double ***f_level = ((FixRespa *) modify->fix[ifix_respa])->f_level;
dtfsq = dtf_inner * step_respa[ilevel];
@ -1298,7 +1312,7 @@ void FixShake::unconstrained_update_respa(int ilevel)
/* ---------------------------------------------------------------------- */
void FixShake::shake2(int m)
void FixShake::shake(int m)
{
int nlist,list[2];
double v[6];

11
src/fix_shake.h
View File

@ -82,6 +82,7 @@ class FixShake : public Fix {
// for angle cluster, 3rd value
// is angletype
double **xshake; // unconstrained atom coords
int *nshake; // count
int vflag; // virial flag
double dtv,dtfsq; // timesteps for trial move
@ -102,13 +103,21 @@ class FixShake : public Fix {
int masscheck(double);
void unconstrained_update();
void unconstrained_update_respa(int);
void shake2(int);
void shake(int);
void shake3(int);
void shake4(int);
void shake3angle(int);
void stats();
int bondfind(int, int, int);
int anglefind(int, int, int);
// static variable for ring communication callback to access class data
// callback functions for ring communication
static FixShake *fsptr;
static void ring_bonds(int, char *);
static void ring_nshake(int, char *);
static void ring_shake(int, char *);
};
}

2
src/read_data.cpp
View File

@ -34,9 +34,9 @@
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "special.h"
#include "error.h"
#include "memory.h"
#include "special.h"
using namespace LAMMPS_NS;

551
src/special.cpp
View File

@ -17,11 +17,16 @@
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
// allocate space for static class variable
Special *Special::sptr;
/* ---------------------------------------------------------------------- */
Special::Special(LAMMPS *lmp) : Pointers(lmp)
@ -49,10 +54,9 @@ Special::~Special()
void Special::build()
{
int i,j,k,m,n,loop,size,original;
int num12,num13,num14;
int max,maxall,messtag,nbuf,nbufmax;
int *buf,*bufcopy,*count;
int i,j,k,m,n,size;
int max,maxall,nbuf;
int *buf;
MPI_Request request;
MPI_Status status;
@ -100,10 +104,7 @@ void Special::build()
nbuf = 0;
for (i = 0; i < nlocal; i++) nbuf += num_bond[i];
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
buf = new int[nbufmax];
bufcopy = new int[nbufmax];
memory->create(buf,nbuf,"special:buf");
// fill buffer with global tags of bond partners of my atoms
@ -116,23 +117,10 @@ void Special::build()
// when receive buffer, scan tags for atoms I own
// when find one, increment nspecial count for that atom
messtag = 1;
for (loop = 0; loop < nprocs; loop++) {
for (i = 0; i < size; i++) {
m = atom->map(buf[i]);
if (m >= 0 && m < nlocal) nspecial[m][0]++;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
sptr = this;
comm->ring(size,sizeof(int),buf,1,ring_one,NULL);
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
}
// ----------------------------------------------------
@ -153,7 +141,7 @@ void Special::build()
// count = accumulating counter
count = new int[nlocal];
memory->create(count,nlocal,"special:count");
for (i = 0; i < nlocal; i++) count[i] = 0;
// add bond partners stored by atom to onetwo list
@ -172,10 +160,7 @@ void Special::build()
nbuf = 0;
for (i = 0; i < nlocal; i++) nbuf += 2*num_bond[i];
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
buf = new int[nbufmax];
bufcopy = new int[nbufmax];
memory->create(buf,nbuf,"special:buf");
// fill buffer with global tags of both atoms in bond
@ -190,26 +175,13 @@ void Special::build()
// when receive buffer, scan 2nd-atom tags for atoms I own
// when find one, add 1st-atom tag to onetwo list for 2nd atom
messtag = 2;
for (loop = 0; loop < nprocs; loop++) {
for (i = 1; i < size; i += 2) {
m = atom->map(buf[i]);
if (m >= 0 && m < nlocal) onetwo[m][count[m]++] = buf[i-1];
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
sptr = this;
comm->ring(size,sizeof(int),buf,2,ring_two,NULL);
memory->destroy(buf);
}
delete [] buf;
delete [] bufcopy;
}
delete [] count;
memory->destroy(count);
// -----------------------------------------------------
// done if special_bonds for 1-3, 1-4 are set to 1.0
@ -230,10 +202,7 @@ void Special::build()
nbuf = 0;
for (i = 0; i < nlocal; i++) nbuf += 2 + nspecial[i][0];
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
buf = new int[nbufmax];
bufcopy = new int[nbufmax];
memory->create(buf,nbuf,"special:buf");
// fill buffer with:
// (1) = counter for 1-3 neighbors, initialized to 0
@ -252,27 +221,8 @@ void Special::build()
// when find one, increment 1-3 count by # of 1-2 neighbors of my atom,
// subtracting one since my list will contain original atom
messtag = 3;
for (loop = 0; loop < nprocs; loop++) {
i = 0;
while (i < size) {
n = buf[i];
num12 = buf[i+1];
for (j = 0; j < num12; j++) {
m = atom->map(buf[i+2+j]);
if (m >= 0 && m < nlocal) n += nspecial[m][0] - 1;
}
buf[i] = n;
i += 2 + num12;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
sptr = this;
comm->ring(size,sizeof(int),buf,3,ring_three,buf);
// extract count from buffer that has cycled back to me
// nspecial[i][1] = # of 1-3 neighbors of atom i
@ -283,8 +233,7 @@ void Special::build()
j += 2 + nspecial[i][0];
}
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
// ----------------------------------------------------
// create onethree[i] = list of 1-3 neighbors for atom i
@ -307,10 +256,7 @@ void Special::build()
nbuf = 0;
for (i = 0; i < nlocal; i++) nbuf += 4 + nspecial[i][0] + nspecial[i][1];
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
buf = new int[nbufmax];
bufcopy = new int[nbufmax];
memory->create(buf,nbuf,"special:buf");
// fill buffer with:
// (1) = global tag of original atom
@ -337,32 +283,8 @@ void Special::build()
// exclude the atom whose tag = original
// this process may include duplicates but they will be culled later
messtag = 4;
for (loop = 0; loop < nprocs; loop++) {
i = 0;
while (i < size) {
original = buf[i];
num12 = buf[i+1];
num13 = buf[i+2];
n = buf[i+3];
for (j = 0; j < num12; j++) {
m = atom->map(buf[i+4+j]);
if (m >= 0 && m < nlocal)
for (k = 0; k < nspecial[m][0]; k++)
if (onetwo[m][k] != original)
buf[i+4+num12+(n++)] = onetwo[m][k];
}
buf[i+3] = n;
i += 4 + num12 + num13;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
sptr = this;
comm->ring(size,sizeof(int),buf,4,ring_four,buf);
// fill onethree with buffer values that have been returned to me
// sanity check: accumulated buf[i+3] count should equal
@ -377,8 +299,7 @@ void Special::build()
onethree[i][k] = buf[j++];
}
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
// done if special_bonds for 1-4 are set to 1.0
@ -397,10 +318,7 @@ void Special::build()
nbuf = 0;
for (i = 0; i < nlocal; i++) nbuf += 2 + nspecial[i][1];
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
buf = new int[nbufmax];
bufcopy = new int[nbufmax];
memory->create(buf,nbuf,"special:buf");
// fill buffer with:
// (1) = counter for 1-4 neighbors, initialized to 0
@ -419,27 +337,8 @@ void Special::build()
// when find one, increment 1-4 count by # of 1-2 neighbors of my atom
// may include duplicates and original atom but they will be culled later
messtag = 5;
for (loop = 0; loop < nprocs; loop++) {
i = 0;
while (i < size) {
n = buf[i];
num13 = buf[i+1];
for (j = 0; j < num13; j++) {
m = atom->map(buf[i+2+j]);
if (m >= 0 && m < nlocal) n += nspecial[m][0];
}
buf[i] = n;
i += 2 + num13;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
sptr = this;
comm->ring(size,sizeof(int),buf,5,ring_five,buf);
// extract count from buffer that has cycled back to me
// nspecial[i][2] = # of 1-4 neighbors of atom i
@ -450,8 +349,7 @@ void Special::build()
j += 2 + nspecial[i][1];
}
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
// ----------------------------------------------------
// create onefour[i] = list of 1-4 neighbors for atom i
@ -475,10 +373,7 @@ void Special::build()
nbuf = 0;
for (i = 0; i < nlocal; i++)
nbuf += 3 + nspecial[i][1] + nspecial[i][2];
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
buf = new int[nbufmax];
bufcopy = new int[nbufmax];
memory->create(buf,nbuf,"special:buf");
// fill buffer with:
// (1) = # of 1-3 neighbors
@ -502,30 +397,8 @@ void Special::build()
// incrementing the count in buf(i+4)
// this process may include duplicates but they will be culled later
messtag = 6;
for (loop = 0; loop < nprocs; loop++) {
i = 0;
while (i < size) {
num13 = buf[i];
num14 = buf[i+1];
n = buf[i+2];
for (j = 0; j < num13; j++) {
m = atom->map(buf[i+3+j]);
if (m >= 0 && m < nlocal)
for (k = 0; k < nspecial[m][0]; k++)
buf[i+3+num13+(n++)] = onetwo[m][k];
}
buf[i+2] = n;
i += 3 + num13 + num14;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
sptr = this;
comm->ring(size,sizeof(int),buf,6,ring_six,buf);
// fill onefour with buffer values that have been returned to me
// sanity check: accumulated buf[i+2] count should equal
@ -540,8 +413,7 @@ void Special::build()
onefour[i][k] = buf[j++];
}
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
combine();
if (force->special_angle) angle_trim();
@ -699,7 +571,7 @@ void Special::combine()
void Special::angle_trim()
{
int i,j,m,n,iglobal,jglobal,ilocal,jlocal;
int i,j,m,n;
MPI_Request request;
MPI_Status status;
@ -743,7 +615,6 @@ void Special::angle_trim()
int maxcount = 0;
for (i = 0; i < nlocal; i++)
maxcount = MAX(maxcount,nspecial[i][1]-nspecial[i][0]);
int **dflag;
memory->create(dflag,nlocal,maxcount,"special::dflag");
for (i = 0; i < nlocal; i++) {
@ -757,12 +628,8 @@ void Special::angle_trim()
int nbuf = 0;
for (i = 0; i < nlocal; i++) nbuf += 2*num_angle[i] + 2*2*num_dihedral[i];
int nbufmax;
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
int *buf = new int[nbufmax];
int *bufcopy = new int[nbufmax];
int *buf;
memory->create(buf,nbuf,"special:buf");
// fill buffer with list of 1,3 atoms in each angle
// and with list of 1,3 and 2,4 atoms in each dihedral
@ -785,41 +652,8 @@ void Special::angle_trim()
// when receive buffer, scan list of 1,3 atoms looking for atoms I own
// when find one, scan its 1-3 neigh list and mark I,J as in an angle
int next = me + 1;
int prev = me -1;
if (next == nprocs) next = 0;
if (prev < 0) prev = nprocs - 1;
int messtag = 7;
for (int loop = 0; loop < nprocs; loop++) {
i = 0;
while (i < size) {
iglobal = buf[i];
jglobal = buf[i+1];
ilocal = atom->map(iglobal);
jlocal = atom->map(jglobal);
if (ilocal >= 0 && ilocal < nlocal)
for (m = nspecial[ilocal][0]; m < nspecial[ilocal][1]; m++)
if (jglobal == special[ilocal][m]) {
dflag[ilocal][m-nspecial[ilocal][0]] = 1;
break;
}
if (jlocal >= 0 && jlocal < nlocal)
for (m = nspecial[jlocal][0]; m < nspecial[jlocal][1]; m++)
if (iglobal == special[jlocal][m]) {
dflag[jlocal][m-nspecial[jlocal][0]] = 1;
break;
}
i += 2;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
sptr = this;
comm->ring(size,sizeof(int),buf,7,ring_seven,NULL);
// delete 1-3 neighbors if they are not flagged in dflag
// preserve 1-4 neighbors
@ -839,8 +673,7 @@ void Special::angle_trim()
// clean up
memory->destroy(dflag);
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
// if no angles or dihedrals are defined,
// delete all 1-3 neighs, preserving 1-4 neighs
@ -880,7 +713,7 @@ void Special::angle_trim()
void Special::dihedral_trim()
{
int i,j,m,n,iglobal,jglobal,ilocal,jlocal;
int i,j,m,n;
MPI_Request request;
MPI_Status status;
@ -918,7 +751,6 @@ void Special::dihedral_trim()
int maxcount = 0;
for (i = 0; i < nlocal; i++)
maxcount = MAX(maxcount,nspecial[i][2]-nspecial[i][1]);
int **dflag;
memory->create(dflag,nlocal,maxcount,"special::dflag");
for (i = 0; i < nlocal; i++) {
@ -931,12 +763,8 @@ void Special::dihedral_trim()
int nbuf = 0;
for (i = 0; i < nlocal; i++) nbuf += 2*num_dihedral[i];
int nbufmax;
MPI_Allreduce(&nbuf,&nbufmax,1,MPI_INT,MPI_MAX,world);
int *buf = new int[nbufmax];
int *bufcopy = new int[nbufmax];
int *buf;
memory->create(buf,nbuf,"special:buf");
// fill buffer with list of 1,4 atoms in each dihedral
@ -951,41 +779,8 @@ void Special::dihedral_trim()
// when receive buffer, scan list of 1,4 atoms looking for atoms I own
// when find one, scan its 1-4 neigh list and mark I,J as in a dihedral
int next = me + 1;
int prev = me -1;
if (next == nprocs) next = 0;
if (prev < 0) prev = nprocs - 1;
int messtag = 8;
for (int loop = 0; loop < nprocs; loop++) {
i = 0;
while (i < size) {
iglobal = buf[i];
jglobal = buf[i+1];
ilocal = atom->map(iglobal);
jlocal = atom->map(jglobal);
if (ilocal >= 0 && ilocal < nlocal)
for (m = nspecial[ilocal][1]; m < nspecial[ilocal][2]; m++)
if (jglobal == special[ilocal][m]) {
dflag[ilocal][m-nspecial[ilocal][1]] = 1;
break;
}
if (jlocal >= 0 && jlocal < nlocal)
for (m = nspecial[jlocal][1]; m < nspecial[jlocal][2]; m++)
if (iglobal == special[jlocal][m]) {
dflag[jlocal][m-nspecial[jlocal][1]] = 1;
break;
}
i += 2;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
MPI_Send(buf,size,MPI_INT,next,messtag,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_INT,&size);
for (j = 0; j < size; j++) buf[j] = bufcopy[j];
}
}
sptr = this;
comm->ring(size,sizeof(int),buf,8,ring_eight,NULL);
// delete 1-4 neighbors if they are not flagged in dflag
@ -1000,8 +795,7 @@ void Special::dihedral_trim()
// clean up
memory->destroy(dflag);
delete [] buf;
delete [] bufcopy;
memory->destroy(buf);
// if no dihedrals are defined, delete all 1-4 neighs
@ -1024,3 +818,254 @@ void Special::dihedral_trim()
" %g = # of 1-4 neighbors after dihedral trim\n",allcount);
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan tags for atoms I own
when find one, increment nspecial count for that atom
------------------------------------------------------------------------- */
void Special::ring_one(int ndatum, char *cbuf)
{
Atom *atom = sptr->atom;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int *buf = (int *) cbuf;
int m;
for (int i = 0; i < ndatum; i++) {
m = atom->map(buf[i]);
if (m >= 0 && m < nlocal) nspecial[m][0]++;
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan 2nd-atom tags for atoms I own
when find one, add 1st-atom tag to onetwo list for 2nd atom
------------------------------------------------------------------------- */
void Special::ring_two(int ndatum, char *cbuf)
{
Atom *atom = sptr->atom;
int nlocal = atom->nlocal;
int **onetwo = sptr->onetwo;
int *count = sptr->count;
int *buf = (int *) cbuf;
int m;
for (int i = 1; i < ndatum; i += 2) {
m = atom->map(buf[i]);
if (m >= 0 && m < nlocal) onetwo[m][count[m]++] = buf[i-1];
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan list of 1-2 neighbors for atoms I own
when find one, increment 1-3 count by # of 1-2 neighbors of my atom,
subtracting one since my list will contain original atom
------------------------------------------------------------------------- */
void Special::ring_three(int ndatum, char *cbuf)
{
Atom *atom = sptr->atom;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int *buf = (int *) cbuf;
int i,j,m,n,num12;
i = 0;
while (i < ndatum) {
n = buf[i];
num12 = buf[i+1];
for (j = 0; j < num12; j++) {
m = atom->map(buf[i+2+j]);
if (m >= 0 && m < nlocal) n += nspecial[m][0] - 1;
}
buf[i] = n;
i += 2 + num12;
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan list of 1-2 neighbors for atoms I own
when find one, add its neighbors to 1-3 list
increment the count in buf(i+4)
exclude the atom whose tag = original
this process may include duplicates but they will be culled later
------------------------------------------------------------------------- */
void Special::ring_four(int ndatum, char *cbuf)
{
Atom *atom = sptr->atom;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int **onetwo = sptr->onetwo;
int *buf = (int *) cbuf;
int i,j,k,m,n,original,num12,num13;
i = 0;
while (i < ndatum) {
original = buf[i];
num12 = buf[i+1];
num13 = buf[i+2];
n = buf[i+3];
for (j = 0; j < num12; j++) {
m = atom->map(buf[i+4+j]);
if (m >= 0 && m < nlocal)
for (k = 0; k < nspecial[m][0]; k++)
if (onetwo[m][k] != original)
buf[i+4+num12+(n++)] = onetwo[m][k];
}
buf[i+3] = n;
i += 4 + num12 + num13;
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan list of 1-3 neighbors for atoms I own
when find one, increment 1-4 count by # of 1-2 neighbors of my atom
may include duplicates and original atom but they will be culled later
------------------------------------------------------------------------- */
void Special::ring_five(int ndatum, char *cbuf)
{
Atom *atom = sptr->atom;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int *buf = (int *) cbuf;
int i,j,m,n,num13;
i = 0;
while (i < ndatum) {
n = buf[i];
num13 = buf[i+1];
for (j = 0; j < num13; j++) {
m = atom->map(buf[i+2+j]);
if (m >= 0 && m < nlocal) n += nspecial[m][0];
}
buf[i] = n;
i += 2 + num13;
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan list of 1-3 neighbors for atoms I own
when find one, add its neighbors to 1-4 list
incrementing the count in buf(i+4)
this process may include duplicates but they will be culled later
------------------------------------------------------------------------- */
void Special::ring_six(int ndatum, char *cbuf)
{
Atom *atom = sptr->atom;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int **onetwo = sptr->onetwo;
int *buf = (int *) cbuf;
int i,j,k,m,n,num13,num14;
i = 0;
while (i < ndatum) {
num13 = buf[i];
num14 = buf[i+1];
n = buf[i+2];
for (j = 0; j < num13; j++) {
m = atom->map(buf[i+3+j]);
if (m >= 0 && m < nlocal)
for (k = 0; k < nspecial[m][0]; k++)
buf[i+3+num13+(n++)] = onetwo[m][k];
}
buf[i+2] = n;
i += 3 + num13 + num14;
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan list of 1,3 atoms looking for atoms I own
when find one, scan its 1-3 neigh list and mark I,J as in an angle
------------------------------------------------------------------------- */
void Special::ring_seven(int ndatum, char *cbuf)
{
Atom *atom = sptr->atom;
int **nspecial = atom->nspecial;
int **special = atom->special;
int nlocal = atom->nlocal;
int **dflag = sptr->dflag;
int *buf = (int *) cbuf;
int i,m,iglobal,jglobal,ilocal,jlocal;
i = 0;
while (i < ndatum) {
iglobal = buf[i];
jglobal = buf[i+1];
ilocal = atom->map(iglobal);
jlocal = atom->map(jglobal);
if (ilocal >= 0 && ilocal < nlocal)
for (m = nspecial[ilocal][0]; m < nspecial[ilocal][1]; m++)
if (jglobal == special[ilocal][m]) {
dflag[ilocal][m-nspecial[ilocal][0]] = 1;
break;
}
if (jlocal >= 0 && jlocal < nlocal)
for (m = nspecial[jlocal][0]; m < nspecial[jlocal][1]; m++)
if (iglobal == special[jlocal][m]) {
dflag[jlocal][m-nspecial[jlocal][0]] = 1;
break;
}
i += 2;
}
}
/* ----------------------------------------------------------------------
when receive buffer, scan list of 1,4 atoms looking for atoms I own
when find one, scan its 1-4 neigh list and mark I,J as in a dihedral
------------------------------------------------------------------------- */
void Special::ring_eight(int ndatum, char *cbuf)
{
Atom *atom = sptr->atom;
int **nspecial = atom->nspecial;
int **special = atom->special;
int nlocal = atom->nlocal;
int **dflag = sptr->dflag;
int *buf = (int *) cbuf;
int i,m,iglobal,jglobal,ilocal,jlocal;
i = 0;
while (i < ndatum) {
iglobal = buf[i];
jglobal = buf[i+1];
ilocal = atom->map(iglobal);
jlocal = atom->map(jglobal);
if (ilocal >= 0 && ilocal < nlocal)
for (m = nspecial[ilocal][1]; m < nspecial[ilocal][2]; m++)
if (jglobal == special[ilocal][m]) {
dflag[ilocal][m-nspecial[ilocal][1]] = 1;
break;
}
if (jlocal >= 0 && jlocal < nlocal)
for (m = nspecial[jlocal][1]; m < nspecial[jlocal][2]; m++)
if (iglobal == special[jlocal][m]) {
dflag[jlocal][m-nspecial[jlocal][1]] = 1;
break;
}
i += 2;
}
}

18
src/special.h
View File

@ -29,9 +29,27 @@ class Special : protected Pointers {
int **onetwo,**onethree,**onefour;
int dihedral_flag;
// data used by ring callback methods
int *count;
int **dflag;
void combine();
void angle_trim();
void dihedral_trim();
// static variable for ring communication callback to access class data
// callback functions for ring communication
static Special *sptr;
static void ring_one(int, char *);
static void ring_two(int, char *);
static void ring_three(int, char *);
static void ring_four(int, char *);
static void ring_five(int, char *);
static void ring_six(int, char *);
static void ring_seven(int, char *);
static void ring_eight(int, char *);
};
}