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@9416 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2013-02-08 21:38:50 +00:00
parent 906a57d529
commit 7a286596d8
12 changed files with 166 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
src/KSPACE/pair_buck_long_coul_long.cpp
View File

@ -89,17 +89,17 @@ void PairBuckLongCoulLong::settings(int narg, char **arg)
ewald_order = 0;
ewald_off = 0;
options(arg, 6);
options(++arg, 1);
options(arg,6);
options(++arg,1);
if (!comm->me && ewald_order&(1<<6)) error->warning(FLERR,PAIR_MIX);
if (!comm->me && ewald_order==((1<<1)|(1<<6)))
if (!comm->me && ewald_order & (1<<6)) error->warning(FLERR,PAIR_MIX);
if (!comm->me && ewald_order == ((1<<1) | (1<<6)))
error->warning(FLERR,PAIR_LARGEST);
if (!*(++arg)) error->all(FLERR,PAIR_MISSING);
if (ewald_off&(1<<6)) error->all(FLERR,PAIR_LJ_OFF);
if (!((ewald_order^ewald_off)&(1<<1))) error->all(FLERR,PAIR_COUL_CUT);
if (ewald_off & (1<<6)) error->all(FLERR,PAIR_LJ_OFF);
if (!((ewald_order^ewald_off) & (1<<1))) error->all(FLERR,PAIR_COUL_CUT);
cut_buck_global = force->numeric(*(arg++));
if (*arg&&(ewald_order&0x42==0x42)) error->all(FLERR,PAIR_CUTOFF);
if (*arg && ((ewald_order & 0x42) == 0x42)) error->all(FLERR,PAIR_CUTOFF);
if (narg == 4) cut_coul = force->numeric(*arg);
else cut_coul = cut_buck_global;

8
src/KSPACE/pair_lj_long_coul_long.cpp
View File

@ -88,13 +88,13 @@ void PairLJLongCoulLong::settings(int narg, char **arg)
ewald_order = 0;
options(arg, 6);
options(++arg, 1);
if (!comm->me && ewald_order&(1<<6)) error->warning(FLERR,PAIR_MIX);
if (!comm->me && ewald_order==((1<<1)|(1<<6)))
if (!comm->me && ewald_order & (1<<6)) error->warning(FLERR,PAIR_MIX);
if (!comm->me && ewald_order == ((1<<1) | (1<<6)))
error->warning(FLERR,PAIR_LARGEST);
if (!*(++arg)) error->all(FLERR,PAIR_MISSING);
if (!((ewald_order^ewald_off)&(1<<1))) error->all(FLERR,PAIR_COUL_CUT);
if (!((ewald_order^ewald_off) & (1<<1))) error->all(FLERR,PAIR_COUL_CUT);
cut_lj_global = force->numeric(*(arg++));
if (*arg&&(ewald_order&0x42==0x42)) error->all(FLERR,PAIR_CUTOFF);
if (*arg && ((ewald_order & 0x42) == 0x42)) error->all(FLERR,PAIR_CUTOFF);
if (narg == 4) cut_coul = force->numeric(*arg);
else cut_coul = cut_lj_global;

50
src/comm.cpp
View File

@ -1454,6 +1454,56 @@ void Comm::reverse_comm_dump(Dump *dump)
}
}
/* ----------------------------------------------------------------------
forward communication of N values in array
------------------------------------------------------------------------- */
void Comm::forward_comm_array(int n, double **array)
{
int i,j,k,m,iswap,last;
double *buf;
MPI_Request request;
MPI_Status status;
// check that buf_send and buf_recv are big enough
for (iswap = 0; iswap < nswap; iswap++) {
// pack buffer
m = 0;
for (i = 0; i < sendnum[iswap]; i++) {
j = sendlist[iswap][i];
for (k = 0; k < n; k++)
buf_send[m++] = array[j][k];
}
// exchange with another proc
// if self, set recv buffer to send buffer
if (sendproc[iswap] != me) {
if (recvnum[iswap])
MPI_Irecv(buf_recv,n*recvnum[iswap],MPI_DOUBLE,recvproc[iswap],0,
world,&request);
if (sendnum[iswap])
MPI_Send(buf_send,n*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0,world);
if (recvnum[iswap]) MPI_Wait(&request,&status);
buf = buf_recv;
} else buf = buf_send;
// unpack buffer
m = 0;
last = firstrecv[iswap] + recvnum[iswap];
for (i = firstrecv[iswap]; i < last; i++)
for (k = 0; k < n; k++)
array[i][k] = buf[m++];
}
}
/* ----------------------------------------------------------------------
reverse communication invoked by a Dump
------------------------------------------------------------------------- */

2
src/comm.h
View File

@ -59,6 +59,8 @@ class Comm : protected Pointers {
virtual void reverse_comm_compute(class Compute *); // reverse from a Compute
virtual void forward_comm_dump(class Dump *); // forward comm from a Dump
virtual void reverse_comm_dump(class Dump *); // reverse comm from a Dump
void forward_comm_array(int, double **); // forward comm of array
void ring(int, int, char *, char *, int,
void (*)(int, char *)); // ring communication

95
src/domain.cpp
View File

@ -529,8 +529,77 @@ void Domain::pbc()
}
/* ----------------------------------------------------------------------
check that no pair of atoms in a bonded interaction
are further apart than half a periodic box length
warn if image flags of any bonded atoms are inconsistent
could be a problem when using replicate or fix rigid
------------------------------------------------------------------------- */
void Domain::image_check()
{
int i,j,k;
// only need to check if system is molecular and some dimension is periodic
// if running verlet/split, don't check on KSpace partition since
// it has no ghost atoms and thus bond partners won't exist
if (!atom->molecular) return;
if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return;
if (strcmp(update->integrate_style,"verlet/split") == 0 &&
universe->iworld != 0) return;
// communicate unwrapped position of owned atoms to ghost atoms
double **unwrap;
memory->create(unwrap,atom->nmax,3,"domain:unwrap");
double **x = atom->x;
int *image = atom->image;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
unmap(x[i],image[i],unwrap[i]);
comm->forward_comm_array(3,unwrap);
// compute unwrapped extent of each bond
// flag if any bond component is longer than 1/2 of periodic box length
// flag if any bond component is longer than non-periodic box length
// which means image flags in that dimension were different
int *num_bond = atom->num_bond;
int **bond_atom = atom->bond_atom;
double delx,dely,delz;
int flag = 0;
for (i = 0; i < nlocal; i++)
for (j = 0; j < num_bond[i]; j++) {
k = atom->map(bond_atom[i][j]);
if (k == -1) error->one(FLERR,"Bond atom missing in image check");
delx = unwrap[i][0] - unwrap[k][0];
dely = unwrap[i][1] - unwrap[k][1];
delz = unwrap[i][2] - unwrap[k][2];
if (xperiodic && delx > xprd_half) flag = 1;
if (xperiodic && dely > yprd_half) flag = 1;
if (dimension == 3 && zperiodic && delz > zprd_half) flag = 1;
if (!xperiodic && delx > xprd) flag = 1;
if (!yperiodic && dely > yprd) flag = 1;
if (dimension == 3 && !zperiodic && delz > zprd) flag = 1;
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flag && comm->me == 0) error->warning(FLERR,"Inconsistent image flags");
memory->destroy(unwrap);
}
/* ----------------------------------------------------------------------
warn if end atoms in any bonded interaction
are further apart than half a periodic box length
could cause problems when bonded neighbor list is built since
closest_image() could return wrong image
------------------------------------------------------------------------- */
void Domain::box_too_small_check()
@ -547,9 +616,10 @@ void Domain::box_too_small_check()
universe->iworld != 0) return;
// maxbondall = longest current bond length
// NOTE: if box is tiny (less than 2 * bond-length),
// the check itself may compute bad bond lengths
// not sure how to account for that extreme case
// if periodic box dim is tiny (less than 2 * bond-length),
// minimum_image() itself may compute bad bond lengths
// in this case, image_check() should warn,
// assuming 2 atoms have consistent image flags
int *num_bond = atom->num_bond;
int **bond_atom = atom->bond_atom;
@ -566,7 +636,7 @@ void Domain::box_too_small_check()
delx = x[i][0] - x[k][0];
dely = x[i][1] - x[k][1];
delz = x[i][2] - x[k][2];
domain->minimum_image(delx,dely,delz);
minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
maxbondme = MAX(maxbondme,rsq);
}
@ -582,18 +652,19 @@ void Domain::box_too_small_check()
if (atom->nangles) maxdelta = 2.0 * maxbondall * BONDSTRETCH;
if (atom->ndihedrals) maxdelta = 3.0 * maxbondall * BONDSTRETCH;
// maxdelta cannot be more than half a periodic box length,
// else when use minimg() in bond/angle/dihdral compute,
// could calculate incorrect distance between 2 atoms
// warn if maxdelta > than half any periodic box length
// since atoms in the interaction could rotate into that dimension
int flag = 0;
if (xperiodic && maxdelta > xprd_half) flag = 1;
if (yperiodic && maxdelta > yprd_half) flag = 1;
if (dimension == 3 && zperiodic && maxdelta > zprd_half) flag = 1;
if (flag)
error->all(FLERR,
"Bond/angle/dihedral extent > half of periodic box length");
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flag && comm->me == 0)
error->warning(FLERR,
"Bond/angle/dihedral extent > half of periodic box length");
}
/* ----------------------------------------------------------------------

1
src/domain.h
View File

@ -95,6 +95,7 @@ class Domain : protected Pointers {
virtual void set_local_box();
virtual void reset_box();
virtual void pbc();
void image_check();
void box_too_small_check();
void minimum_image(double &, double &, double &);
void minimum_image(double *);

1
src/fix.h
View File

@ -125,6 +125,7 @@ class Fix : protected Pointers {
virtual void min_setup_pre_neighbor() {}
virtual void min_setup_pre_force(int) {}
virtual void min_pre_exchange() {}
virtual void min_pre_neighbor() {}
virtual void min_pre_force(int) {}
virtual void min_post_force(int) {}

3
src/min.cpp
View File

@ -234,6 +234,7 @@ void Min::setup()
if (atom->sortfreq > 0) atom->sort();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
modify->setup_pre_neighbor();
neighbor->build();
@ -321,6 +322,7 @@ void Min::setup_minimal(int flag)
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
modify->setup_pre_neighbor();
neighbor->build();
@ -436,7 +438,6 @@ void Min::cleanup()
// delete fix at end of run, so its atom arrays won't persist
modify->delete_fix("MINIMIZE");
domain->box_too_small_check();
}

12
src/modify.cpp
View File

@ -56,7 +56,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
list_initial_integrate_respa = list_post_integrate_respa = NULL;
list_pre_force_respa = list_post_force_respa = NULL;
list_final_integrate_respa = NULL;
list_min_pre_exchange = list_pre_neighbor = NULL;
list_min_pre_exchange = list_min_pre_neighbor = NULL;
list_min_pre_force = list_min_post_force = NULL;
list_min_energy = NULL;
@ -454,6 +454,16 @@ void Modify::min_pre_exchange()
fix[list_min_pre_exchange[i]]->min_pre_exchange();
}
/* ----------------------------------------------------------------------
minimizer pre-neighbor call, only for relevant fixes
------------------------------------------------------------------------- */
void Modify::min_pre_neighbor()
{
for (int i = 0; i < n_min_pre_neighbor; i++)
fix[list_min_pre_neighbor[i]]->min_pre_neighbor();
}
/* ----------------------------------------------------------------------
minimizer pre-force call, only for relevant fixes
------------------------------------------------------------------------- */

1
src/modify.h
View File

@ -67,6 +67,7 @@ class Modify : protected Pointers {
void final_integrate_respa(int, int);
void min_pre_exchange();
void min_pre_neighbor();
void min_pre_force(int);
void min_post_force(int);

2
src/respa.cpp
View File

@ -351,6 +351,7 @@ void Respa::setup()
if (atom->sortfreq > 0) atom->sort();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
modify->setup_pre_neighbor();
neighbor->build();
@ -417,6 +418,7 @@ void Respa::setup_minimal(int flag)
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
modify->setup_pre_neighbor();
neighbor->build();

2
src/verlet.cpp
View File

@ -111,6 +111,7 @@ void Verlet::setup()
if (atom->sortfreq > 0) atom->sort();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
modify->setup_pre_neighbor();
neighbor->build();
@ -169,6 +170,7 @@ void Verlet::setup_minimal(int flag)
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
modify->setup_pre_neighbor();
neighbor->build();