jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2007-03-08 00:54:02 +00:00
parent a70ee2eadc
commit 9e1f8bcd6a
64 changed files with 1785 additions and 772 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
src/atom.cpp
View File

@ -35,18 +35,17 @@
using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
#define DELTA 1
#define DELTA_MEMSTR 1024
#define EPSILON 1.0e-6
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
/* ---------------------------------------------------------------------- */
Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
{
// no atoms
natoms = nlocal = nghost = nmax = 0;
ntypes = 0;
nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
@ -620,8 +619,9 @@ int Atom::count_words(char *line)
void Atom::data_atoms(int n, char *buf, int ihybrid)
{
int m,imagetmp,xptr,iptr;
double xtmp,ytmp,ztmp;
int m,imagedata,xptr,iptr;
double xdata[3],lamda[3],sublo[3],subhi[3];
double *coord;
char *next;
next = strchr(buf,'\n');
@ -635,12 +635,33 @@ void Atom::data_atoms(int n, char *buf, int ihybrid)
char **values = new char*[nwords];
double subxlo = domain->subxlo;
double subxhi = domain->subxhi;
double subylo = domain->subylo;
double subyhi = domain->subyhi;
double subzlo = domain->subzlo;
double subzhi = domain->subzhi;
// set bounds for my proc
// if periodic and I am lo/hi proc, adjust bounds by EPSILON
// insures all data atoms will be owned even with round-off
int triclinic = domain->triclinic;
if (triclinic == 0) {
sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0];
sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1];
sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2];
} else {
sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0];
sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1];
sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2];
}
if (domain->xperiodic) {
if (comm->myloc[0] == 0) sublo[0] -= EPSILON;
if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += EPSILON;
}
if (domain->yperiodic) {
if (comm->myloc[1] == 0) sublo[1] -= EPSILON;
if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += EPSILON;
}
if (domain->zperiodic) {
if (comm->myloc[2] == 0) sublo[2] -= EPSILON;
if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += EPSILON;
}
// xptr = which word in line starts xyz coords
// iptr = which word in line starts ix,iy,iz image flags
@ -662,20 +683,25 @@ void Atom::data_atoms(int n, char *buf, int ihybrid)
for (m = 1; m < nwords; m++)
values[m] = strtok(NULL," \t\n\r\f");
xtmp = atof(values[xptr]);
ytmp = atof(values[xptr+1]);
ztmp = atof(values[xptr+2]);
if (imageflag)
imagetmp = ((atoi(values[iptr+2]) + 512 & 1023) << 20) |
imagedata = ((atoi(values[iptr+2]) + 512 & 1023) << 20) |
((atoi(values[iptr+1]) + 512 & 1023) << 10) |
(atoi(values[iptr]) + 512 & 1023);
else imagetmp = (512 << 20) | (512 << 10) | 512;
else imagedata = (512 << 20) | (512 << 10) | 512;
domain->remap(xtmp,ytmp,ztmp,imagetmp);
if (xtmp >= subxlo && xtmp < subxhi &&
ytmp >= subylo && ytmp < subyhi &&
ztmp >= subzlo && ztmp < subzhi)
avec->data_atom(xtmp,ytmp,ztmp,imagetmp,values,ihybrid);
xdata[0] = atof(values[xptr]);
xdata[1] = atof(values[xptr+1]);
xdata[2] = atof(values[xptr+2]);
domain->remap(xdata,imagedata);
if (triclinic) {
domain->x2lamda(xdata,lamda);
coord = lamda;
} else coord = xdata;
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2])
avec->data_atom(xdata,imagedata,values,ihybrid);
buf = next + 1;
}

8
src/atom_vec.h
View File

@ -44,7 +44,7 @@ class AtomVec : protected Pointers {
virtual void zero_ghost(int,int) {}
virtual void copy(int, int) = 0;
virtual int pack_comm(int, int *, double *, int *) = 0;
virtual int pack_comm(int, int *, double *, int, double *) = 0;
virtual int pack_comm_one(int, double *) {return 0;}
virtual void unpack_comm(int, int, double *) = 0;
virtual int unpack_comm_one(int, double *) {return 0;}
@ -54,7 +54,7 @@ class AtomVec : protected Pointers {
virtual void unpack_reverse(int, int *, double *) = 0;
virtual int unpack_reverse_one(int, double *) {return 0;}
virtual int pack_border(int, int *, double *, int *) = 0;
virtual int pack_border(int, int *, double *, int, double *) = 0;
virtual int pack_border_one(int, double *) {return 0;}
virtual void unpack_border(int, int, double *) = 0;
virtual int unpack_border_one(int, double *) {return 0;}
@ -67,8 +67,8 @@ class AtomVec : protected Pointers {
virtual int pack_restart(int, double *) = 0;
virtual int unpack_restart(double *) = 0;
virtual void create_atom(int, double, double, double, int) = 0;
virtual void data_atom(double, double, double, int, char **, int) = 0;
virtual void create_atom(int, double *, int) = 0;
virtual void data_atom(double *, int, char **, int) = 0;
virtual void data_vel(int, char *, int);
virtual void data_params(int) {}

50
src/atom_vec_atomic.cpp
View File

@ -14,7 +14,6 @@
#include "stdlib.h"
#include "atom_vec_atomic.h"
#include "atom.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "memory.h"
@ -102,12 +101,13 @@ void AtomVecAtomic::copy(int i, int j)
/* ---------------------------------------------------------------------- */
int AtomVecAtomic::pack_comm(int n, int *list, double *buf, int *pbc_flags)
int AtomVecAtomic::pack_comm(int n, int *list, double *buf,
int pbc_flag, double *pbc_dist)
{
int i,j,m;
m = 0;
if (pbc_flags[0] == 0) {
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
@ -115,14 +115,11 @@ int AtomVecAtomic::pack_comm(int n, int *list, double *buf, int *pbc_flags)
buf[m++] = x[j][2];
}
} else {
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + pbc_flags[1]*xprd;
buf[m++] = x[j][1] + pbc_flags[2]*yprd;
buf[m++] = x[j][2] + pbc_flags[3]*zprd;
buf[m++] = x[j][0] + pbc_dist[0];
buf[m++] = x[j][1] + pbc_dist[1];
buf[m++] = x[j][2] + pbc_dist[2];
}
}
return m;
@ -176,12 +173,13 @@ void AtomVecAtomic::unpack_reverse(int n, int *list, double *buf)
/* ---------------------------------------------------------------------- */
int AtomVecAtomic::pack_border(int n, int *list, double *buf, int *pbc_flags)
int AtomVecAtomic::pack_border(int n, int *list, double *buf,
int pbc_flag, double *pbc_dist)
{
int i,j,m;
m = 0;
if (pbc_flags[0] == 0) {
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
@ -192,14 +190,11 @@ int AtomVecAtomic::pack_border(int n, int *list, double *buf, int *pbc_flags)
buf[m++] = mask[j];
}
} else {
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + pbc_flags[1]*xprd;
buf[m++] = x[j][1] + pbc_flags[2]*yprd;
buf[m++] = x[j][2] + pbc_flags[3]*zprd;
buf[m++] = x[j][0] + pbc_dist[0];
buf[m++] = x[j][1] + pbc_dist[1];
buf[m++] = x[j][2] + pbc_dist[2];
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
@ -378,21 +373,20 @@ int AtomVecAtomic::unpack_restart(double *buf)
}
/* ----------------------------------------------------------------------
create one atom of type itype at x0,y0,z0
create one atom of itype at coord
set other values to defaults
------------------------------------------------------------------------- */
void AtomVecAtomic::create_atom(int itype, double x0, double y0, double z0,
int ihybrid)
void AtomVecAtomic::create_atom(int itype, double *coord, int ihybrid)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
tag[nlocal] = 0;
type[nlocal] = itype;
x[nlocal][0] = x0;
x[nlocal][1] = y0;
x[nlocal][2] = z0;
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = (512 << 20) | (512 << 10) | 512;
v[nlocal][0] = 0.0;
@ -407,8 +401,8 @@ void AtomVecAtomic::create_atom(int itype, double x0, double y0, double z0,
initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecAtomic::data_atom(double xtmp, double ytmp, double ztmp,
int imagetmp, char **values, int ihybrid)
void AtomVecAtomic::data_atom(double *coord, int imagetmp, char **values,
int ihybrid)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
@ -421,9 +415,9 @@ void AtomVecAtomic::data_atom(double xtmp, double ytmp, double ztmp,
if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
error->one("Invalid atom type in Atoms section of data file");
x[nlocal][0] = xtmp;
x[nlocal][1] = ytmp;
x[nlocal][2] = ztmp;
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
image[nlocal] = imagetmp;

8
src/atom_vec_atomic.h
View File

@ -25,11 +25,11 @@ class AtomVecAtomic : public AtomVec {
void grow(int);
void reset_ptrs();
void copy(int, int);
virtual int pack_comm(int, int *, double *, int *);
virtual int pack_comm(int, int *, double *, int, double *);
virtual void unpack_comm(int, int, double *);
int pack_reverse(int, int, double *);
void unpack_reverse(int, int *, double *);
virtual int pack_border(int, int *, double *, int *);
virtual int pack_border(int, int *, double *, int, double *);
virtual void unpack_border(int, int, double *);
int pack_exchange(int, double *);
int unpack_exchange(double *);
@ -37,8 +37,8 @@ class AtomVecAtomic : public AtomVec {
int size_restart_one(int);
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double, double, double, int);
void data_atom(double, double, double, int, char **, int);
void create_atom(int, double *, int);
void data_atom(double *, int, char **, int);
int memory_usage();
protected:

50
src/atom_vec_charge.cpp
View File

@ -14,7 +14,6 @@
#include "stdlib.h"
#include "atom_vec_charge.h"
#include "atom.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "memory.h"
@ -134,12 +133,13 @@ void AtomVecCharge::copy(int i, int j)
/* ---------------------------------------------------------------------- */
int AtomVecCharge::pack_comm(int n, int *list, double *buf, int *pbc_flags)
int AtomVecCharge::pack_comm(int n, int *list, double *buf,
int pbc_flag, double *pbc_dist)
{
int i,j,m;
m = 0;
if (pbc_flags[0] == 0) {
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
@ -147,14 +147,11 @@ int AtomVecCharge::pack_comm(int n, int *list, double *buf, int *pbc_flags)
buf[m++] = x[j][2];
}
} else {
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + pbc_flags[1]*xprd;
buf[m++] = x[j][1] + pbc_flags[2]*yprd;
buf[m++] = x[j][2] + pbc_flags[3]*zprd;
buf[m++] = x[j][0] + pbc_dist[0];
buf[m++] = x[j][1] + pbc_dist[1];
buf[m++] = x[j][2] + pbc_dist[2];
}
}
return m;
@ -208,12 +205,13 @@ void AtomVecCharge::unpack_reverse(int n, int *list, double *buf)
/* ---------------------------------------------------------------------- */
int AtomVecCharge::pack_border(int n, int *list, double *buf, int *pbc_flags)
int AtomVecCharge::pack_border(int n, int *list, double *buf,
int pbc_flag, double *pbc_dist)
{
int i,j,m;
m = 0;
if (pbc_flags[0] == 0) {
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
@ -225,14 +223,11 @@ int AtomVecCharge::pack_border(int n, int *list, double *buf, int *pbc_flags)
buf[m++] = q[j];
}
} else {
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + pbc_flags[1]*xprd;
buf[m++] = x[j][1] + pbc_flags[2]*yprd;
buf[m++] = x[j][2] + pbc_flags[3]*zprd;
buf[m++] = x[j][0] + pbc_dist[0];
buf[m++] = x[j][1] + pbc_dist[1];
buf[m++] = x[j][2] + pbc_dist[2];
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
@ -437,21 +432,20 @@ int AtomVecCharge::unpack_restart(double *buf)
}
/* ----------------------------------------------------------------------
create one atom of type itype at x0,y0,z0
create one atom of itype at coord
set other values to defaults
------------------------------------------------------------------------- */
void AtomVecCharge::create_atom(int itype, double x0, double y0, double z0,
int ihybrid)
void AtomVecCharge::create_atom(int itype, double *coord, int ihybrid)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
tag[nlocal] = 0;
type[nlocal] = itype;
x[nlocal][0] = x0;
x[nlocal][1] = y0;
x[nlocal][2] = z0;
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = (512 << 20) | (512 << 10) | 512;
v[nlocal][0] = 0.0;
@ -468,8 +462,8 @@ void AtomVecCharge::create_atom(int itype, double x0, double y0, double z0,
initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecCharge::data_atom(double xtmp, double ytmp, double ztmp,
int imagetmp, char **values, int ihybrid)
void AtomVecCharge::data_atom(double *coord, int imagetmp, char **values,
int ihybrid)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
@ -484,9 +478,9 @@ void AtomVecCharge::data_atom(double xtmp, double ytmp, double ztmp,
q[nlocal] = atof(values[2]);
x[nlocal][0] = xtmp;
x[nlocal][1] = ytmp;
x[nlocal][2] = ztmp;
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
image[nlocal] = imagetmp;

8
src/atom_vec_charge.h
View File

@ -26,11 +26,11 @@ class AtomVecCharge : public AtomVec {
void zero_owned(int);
void zero_ghost(int, int);
void copy(int, int);
int pack_comm(int, int *, double *, int *);
int pack_comm(int, int *, double *, int, double *);
void unpack_comm(int, int, double *);
int pack_reverse(int, int, double *);
void unpack_reverse(int, int *, double *);
int pack_border(int, int *, double *, int *);
int pack_border(int, int *, double *, int, double *);
int pack_border_one(int, double *);
void unpack_border(int, int, double *);
int unpack_border_one(int, double *);
@ -40,8 +40,8 @@ class AtomVecCharge : public AtomVec {
int size_restart_one(int);
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double, double, double, int);
void data_atom(double, double, double, int, char **, int);
void create_atom(int, double *, int);
void data_atom(double *, int, char **, int);
int memory_usage();
private:

43
src/atom_vec_hybrid.cpp
View File

@ -1,3 +1,4 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
@ -14,7 +15,6 @@
#include "string.h"
#include "atom_vec_hybrid.h"
#include "atom.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "memory.h"
@ -152,12 +152,13 @@ void AtomVecHybrid::copy(int i, int j)
/* ---------------------------------------------------------------------- */
int AtomVecHybrid::pack_comm(int n, int *list, double *buf, int *pbc_flags)
int AtomVecHybrid::pack_comm(int n, int *list, double *buf,
int pbc_flag, double *pbc_dist)
{
int i,j,m;
m = 0;
if (pbc_flags[0] == 0) {
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
@ -166,14 +167,11 @@ int AtomVecHybrid::pack_comm(int n, int *list, double *buf, int *pbc_flags)
m += styles[hybrid[j]]->pack_comm_one(j,&buf[m]);
}
} else {
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + pbc_flags[1]*xprd;
buf[m++] = x[j][1] + pbc_flags[2]*yprd;
buf[m++] = x[j][2] + pbc_flags[3]*zprd;
buf[m++] = x[j][0] + pbc_dist[0];
buf[m++] = x[j][1] + pbc_dist[1];
buf[m++] = x[j][2] + pbc_dist[2];
m += styles[hybrid[j]]->pack_comm_one(j,&buf[m]);
}
}
@ -231,12 +229,13 @@ void AtomVecHybrid::unpack_reverse(int n, int *list, double *buf)
/* ---------------------------------------------------------------------- */
int AtomVecHybrid::pack_border(int n, int *list, double *buf, int *pbc_flags)
int AtomVecHybrid::pack_border(int n, int *list, double *buf,
int pbc_flag, double *pbc_dist)
{
int i,j,m;
m = 0;
if (pbc_flags[0] == 0) {
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
@ -249,14 +248,11 @@ int AtomVecHybrid::pack_border(int n, int *list, double *buf, int *pbc_flags)
m += styles[hybrid[j]]->pack_border_one(j,&buf[m]);
}
} else {
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + pbc_flags[1]*xprd;
buf[m++] = x[j][1] + pbc_flags[2]*yprd;
buf[m++] = x[j][2] + pbc_flags[3]*zprd;
buf[m++] = x[j][0] + pbc_dist[0];
buf[m++] = x[j][1] + pbc_dist[1];
buf[m++] = x[j][2] + pbc_dist[2];
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
@ -396,14 +392,13 @@ int AtomVecHybrid::unpack_restart(double *buf)
}
/* ----------------------------------------------------------------------
create one atom of type itype at x0,y0,z0 for ihybrid style
create one atom of itype at coord for ihybrid style
sub-style does create
grow() must occur here so arrays for all sub-styles are grown
zero auxiliary arrays for all other styles before create
------------------------------------------------------------------------- */
void AtomVecHybrid::create_atom(int itype, double x0, double y0, double z0,
int ihybrid)
void AtomVecHybrid::create_atom(int itype, double *coord, int ihybrid)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
@ -411,7 +406,7 @@ void AtomVecHybrid::create_atom(int itype, double x0, double y0, double z0,
for (int m = 0; m < nstyles; m++)
if (m != ihybrid) styles[m]->zero_owned(nlocal);
hybrid[nlocal] = ihybrid;
styles[ihybrid]->create_atom(itype,x0,y0,z0,0);
styles[ihybrid]->create_atom(itype,coord,0);
}
/* ----------------------------------------------------------------------
@ -420,8 +415,8 @@ void AtomVecHybrid::create_atom(int itype, double x0, double y0, double z0,
sub-style will increment nlocal
------------------------------------------------------------------------- */
void AtomVecHybrid::data_atom(double xtmp, double ytmp, double ztmp,
int imagetmp, char **values, int ihybrid)
void AtomVecHybrid::data_atom(double *coord, int imagetmp, char **values,
int ihybrid)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
@ -429,7 +424,7 @@ void AtomVecHybrid::data_atom(double xtmp, double ytmp, double ztmp,
for (int m = 0; m < nstyles; m++)
if (m != ihybrid) styles[m]->zero_owned(nlocal);
hybrid[nlocal] = ihybrid;
styles[ihybrid]->data_atom(xtmp,ytmp,ztmp,imagetmp,values,0);
styles[ihybrid]->data_atom(coord,imagetmp,values,0);
}
/* ----------------------------------------------------------------------

8
src/atom_vec_hybrid.h
View File

@ -29,11 +29,11 @@ class AtomVecHybrid : public AtomVec {
void grow(int);
void reset_ptrs();
void copy(int, int);
int pack_comm(int, int *, double *, int *);
int pack_comm(int, int *, double *, int, double *);
void unpack_comm(int, int, double *);
int pack_reverse(int, int, double *);
void unpack_reverse(int, int *, double *);
int pack_border(int, int *, double *, int *);
int pack_border(int, int *, double *, int, double *);
void unpack_border(int, int, double *);
int pack_exchange(int, double *);
int unpack_exchange(double *);
@ -41,8 +41,8 @@ class AtomVecHybrid : public AtomVec {
int size_restart_one(int) {return 0;}
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double, double, double, int);
void data_atom(double, double, double, int, char **, int);
void create_atom(int, double *, int);
void data_atom(double *, int, char **, int);
void data_vel(int, char *, int);
void data_params(int);
int memory_usage();

200
src/comm.cpp
View File

@ -11,10 +11,11 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "mpi.h"
#include "math.h"
#include "string.h"
#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#include "comm.h"
#include "atom.h"
#include "atom_vec.h"
@ -72,9 +73,7 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp)
maxforward_pair = maxreverse_pair = 0;
}
/* ----------------------------------------------------------------------
destructor, free all memory
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
Comm::~Comm()
{
@ -122,11 +121,13 @@ void Comm::init()
}
/* ----------------------------------------------------------------------
setup 3d grid of procs based on box size
setup 3d grid of procs based on box size
------------------------------------------------------------------------- */
void Comm::set_procs()
{
triclinic = domain->triclinic;
if (user_procgrid[0] == 0) procs2box();
else {
procgrid[0] = user_procgrid[0];
@ -151,6 +152,10 @@ void Comm::set_procs()
MPI_Cart_shift(cartesian,2,1,&procneigh[2][0],&procneigh[2][1]);
MPI_Comm_free(&cartesian);
// can't set lamda box params until procs are assigned
if (triclinic) domain->set_lamda_box();
if (me == 0) {
if (screen) fprintf(screen," %d by %d by %d processor grid\n",
procgrid[0],procgrid[1],procgrid[2]);
@ -166,25 +171,47 @@ void Comm::set_procs()
void Comm::setup()
{
// for triclinic, tilted planes are neigh cutoff apart
// but cutneigh is distance between those planes along xyz (non-tilted) axes
double cutneigh[3];
double *prd,*prd_border,*sublo,*subhi;
if (triclinic == 0) {
prd = prd_border = domain->prd;
sublo = domain->sublo;
subhi = domain->subhi;
cutneigh[0] = cutneigh[1] = cutneigh[2] = neighbor->cutneigh;
} else {
prd = domain->prd;
prd_border = domain->prd_lamda;
sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda;
double *h_inv = domain->h_inv;
double length;
length = sqrt(h_inv[0]*h_inv[0] + h_inv[5]*h_inv[5] + h_inv[4]*h_inv[4]);
cutneigh[0] = neighbor->cutneigh * length;
length = sqrt(h_inv[1]*h_inv[1] + h_inv[3]*h_inv[3]);
cutneigh[1] = neighbor->cutneigh * length;
length = h_inv[2];
cutneigh[2] = neighbor->cutneigh * length;
}
// need = # of procs I need atoms from in each dim
need[0] = static_cast<int>
(neighbor->cutneigh * procgrid[0] / domain->xprd) + 1;
need[1] = static_cast<int>
(neighbor->cutneigh * procgrid[1] / domain->yprd) + 1;
need[2] = static_cast<int>
(neighbor->cutneigh * procgrid[2] / domain->zprd) + 1;
// for 2d, don't communicate in z
need[0] = static_cast<int> (cutneigh[0] * procgrid[0] / prd_border[0]) + 1;
need[1] = static_cast<int> (cutneigh[1] * procgrid[1] / prd_border[1]) + 1;
need[2] = static_cast<int> (cutneigh[2] * procgrid[2] / prd_border[2]) + 1;
if (force->dimension == 2) need[2] = 0;
// if non-periodic, do not communicate further than procgrid-1 away
// this enables very large cutoffs in non-periodic systems
if (domain->xperiodic == 0) need[0] = MIN(need[0],procgrid[0]-1);
if (domain->yperiodic == 0) need[1] = MIN(need[1],procgrid[1]-1);
if (domain->zperiodic == 0) need[2] = MIN(need[2],procgrid[2]-1);
int *periodicity = domain->periodicity;
if (periodicity[0] == 0) need[0] = MIN(need[0],procgrid[0]-1);
if (periodicity[1] == 0) need[1] = MIN(need[1],procgrid[1]-1);
if (periodicity[2] == 0) need[2] = MIN(need[2],procgrid[2]-1);
// allocate comm memory
@ -201,47 +228,62 @@ void Comm::setup()
// set slablo > slabhi for swaps across non-periodic boundaries
// this insures no atoms are swapped
// only for procs owning sub-box at non-periodic end of global box
// pbc_flags = add-on factor for atoms sent across a periodic global boundary
// 0 = not across a boundary
// 1 = add box-length to coord when sending
// -1 = subtract box-length from coord when sending
// pbc_flag: 0 = not across a boundary, 1 = yes across a boundary
// pbc_dist/border: factors to add to atom coords across PBC for comm/borders
// for triclinic, slablo/hi and pbc_border will be used in lamda (0-1) coords
// 1st part of if statement is sending to the west/south/down
// 2nd part of if statement is sending to the east/north/up
int iswap = 0;
for (int dim = 0; dim < 3; dim++) {
for (int ineed = 0; ineed < 2*need[dim]; ineed++) {
pbc_flags[iswap][0] = 0;
pbc_flags[iswap][1] = 0;
pbc_flags[iswap][2] = 0;
pbc_flags[iswap][3] = 0;
pbc_flag[iswap] = 0;
pbc_dist[iswap][0] = pbc_dist[iswap][1] = pbc_dist[iswap][2] = 0.0;
pbc_dist_border[iswap][0] = pbc_dist_border[iswap][1] =
pbc_dist_border[iswap][2] = 0.0;
if (ineed % 2 == 0) {
sendproc[iswap] = procneigh[dim][0];
recvproc[iswap] = procneigh[dim][1];
if (ineed < 2) slablo[iswap] = -BIG;
else slablo[iswap] = 0.5 * (domain->sublo[dim] + domain->subhi[dim]);
slabhi[iswap] = domain->sublo[dim] + neighbor->cutneigh;
else slablo[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
slabhi[iswap] = sublo[dim] + cutneigh[dim];
if (myloc[dim] == 0) {
if (domain->periodicity[dim] == 0)
if (periodicity[dim] == 0)
slabhi[iswap] = slablo[iswap] - 1.0;
else {
pbc_flags[iswap][0] = 1;
pbc_flags[iswap][dim+1] = 1;
pbc_flag[iswap] = 1;
pbc_dist[iswap][dim] = prd[dim];
pbc_dist_border[iswap][dim] = prd_border[dim];
if (triclinic) {
if (dim == 1) pbc_dist[iswap][0] += domain->xy;
else if (dim == 2) {
pbc_dist[iswap][0] += domain->xz;
pbc_dist[iswap][1] += domain->yz;
}
}
}
}
} else {
sendproc[iswap] = procneigh[dim][1];
recvproc[iswap] = procneigh[dim][0];
slablo[iswap] = domain->subhi[dim] - neighbor->cutneigh;
slablo[iswap] = subhi[dim] - cutneigh[dim];
if (ineed < 2) slabhi[iswap] = BIG;
else slabhi[iswap] = 0.5 * (domain->sublo[dim] + domain->subhi[dim]);
else slabhi[iswap] = 0.5 * (sublo[dim] + subhi[dim]);
if (myloc[dim] == procgrid[dim]-1) {
if (domain->periodicity[dim] == 0)
if (periodicity[dim] == 0)
slabhi[iswap] = slablo[iswap] - 1.0;
else {
pbc_flags[iswap][0] = 1;
pbc_flags[iswap][dim+1] = -1;
pbc_flag[iswap] = 1;
pbc_dist[iswap][dim] = -prd[dim];
pbc_dist_border[iswap][dim] = -prd_border[dim];
if (triclinic) {
if (dim == 1) pbc_dist[iswap][0] -= domain->xy;
else if (dim == 2) {
pbc_dist[iswap][0] -= domain->xz;
pbc_dist[iswap][1] -= domain->yz;
}
}
}
}
}
@ -277,14 +319,14 @@ void Comm::communicate()
MPI_Irecv(buf,size_comm_recv[iswap],MPI_DOUBLE,
recvproc[iswap],0,world,&request);
n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
buf_send,pbc_flags[iswap]);
buf_send,pbc_flag[iswap],pbc_dist[iswap]);
MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
MPI_Wait(&request,&status);
} else {
MPI_Irecv(buf_recv,size_comm_recv[iswap],MPI_DOUBLE,
recvproc[iswap],0,world,&request);
n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
buf_send,pbc_flags[iswap]);
buf_send,pbc_flag[iswap],pbc_dist[iswap]);
MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world);
MPI_Wait(&request,&status);
avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_recv);
@ -294,10 +336,11 @@ void Comm::communicate()
if (comm_x_only) {
if (sendnum[iswap])
n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
x[firstrecv[iswap]],pbc_flags[iswap]);
x[firstrecv[iswap]],
pbc_flag[iswap],pbc_dist[iswap]);
} else {
n = avec->pack_comm(sendnum[iswap],sendlist[iswap],
buf_send,pbc_flags[iswap]);
buf_send,pbc_flag[iswap],pbc_dist[iswap]);
avec->unpack_comm(recvnum[iswap],firstrecv[iswap],buf_send);
}
}
@ -363,6 +406,7 @@ void Comm::reverse_communicate()
can happen if atom moves outside of non-periodic bounary
or if atom moves more than one proc away
this routine called before every reneighboring
for triclinic, atoms must be in lamda coords (0-1) before exchange is called
------------------------------------------------------------------------- */
void Comm::exchange()
@ -370,7 +414,7 @@ void Comm::exchange()
int i,m,nsend,nrecv,nrecv1,nrecv2,nlocal;
double lo,hi,value;
double **x;
double *buf;
double *sublo,*subhi,*buf;
MPI_Request request;
MPI_Status status;
AtomVec *avec = atom->avec;
@ -379,6 +423,16 @@ void Comm::exchange()
if (map_style) atom->map_clear();
// subbox bounds for orthogonal or triclinic
if (triclinic == 0) {
sublo = domain->sublo;
subhi = domain->subhi;
} else {
sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda;
}
// loop over dimensions
for (int dim = 0; dim < 3; dim++) {
@ -387,8 +441,8 @@ void Comm::exchange()
// when atom is deleted, fill it in with last atom
x = atom->x;
lo = domain->sublo[dim];
hi = domain->subhi[dim];
lo = sublo[dim];
hi = subhi[dim];
nlocal = atom->nlocal;
i = nsend = 0;
@ -457,6 +511,7 @@ void Comm::exchange()
this does equivalent of a communicate (so don't need to explicitly
call communicate routine on reneighboring timestep)
this routine is called before every reneighboring
for triclinic, atoms must be in lamda coords (0-1) before borders is called
------------------------------------------------------------------------- */
void Comm::borders()
@ -509,7 +564,8 @@ void Comm::borders()
if (nsend*size_border > maxsend)
grow_send(nsend*size_border,0);
n = avec->pack_border(nsend,sendlist[iswap],buf_send,pbc_flags[iswap]);
n = avec->pack_border(nsend,sendlist[iswap],buf_send,
pbc_flag[iswap],pbc_dist_border[iswap]);
// swap atoms with other proc
// put incoming ghosts at end of my atom arrays
@ -577,7 +633,7 @@ void Comm::comm_pair(Pair *pair)
// pack buffer
n = pair->pack_comm(sendnum[iswap],sendlist[iswap],
buf_send,pbc_flags[iswap]);
buf_send,pbc_flag[iswap],pbc_dist[iswap]);
// exchange with another proc
// if self, set recv buffer to send buffer
@ -646,7 +702,7 @@ void Comm::comm_fix(Fix *fix)
// pack buffer
n = fix->pack_comm(sendnum[iswap],sendlist[iswap],
buf_send,pbc_flags[iswap]);
buf_send,pbc_flag[iswap],pbc_dist[iswap]);
// exchange with another proc
// if self, set recv buffer to send buffer
@ -715,7 +771,7 @@ void Comm::comm_compute(Compute *compute)
// pack buffer
n = compute->pack_comm(sendnum[iswap],sendlist[iswap],
buf_send,pbc_flags[iswap]);
buf_send,pbc_flag[iswap],pbc_dist[iswap]);
// exchange with another proc
// if self, set recv buffer to send buffer
@ -770,23 +826,38 @@ void Comm::reverse_comm_compute(Compute *compute)
/* ----------------------------------------------------------------------
assign nprocs to 3d xprd,yprd,zprd box so as to minimize surface area
area = surface area of each of 3 faces of simulation box
for triclinic, area = cross product of 2 edge vectors stored in h matrix
------------------------------------------------------------------------- */
void Comm::procs2box()
{
int ipx,ipy,ipz,nremain;
double boxx,boxy,boxz,surf;
double area[3];
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
if (triclinic == 0) {
area[0] = domain->xprd * domain->yprd;
area[1] = domain->xprd * domain->zprd;
area[2] = domain->yprd * domain->zprd;
} else {
double *h = domain->h;
double x,y,z;
cross(h[0],0.0,0.0,h[5],h[1],0.0,x,y,z);
area[0] = sqrt(x*x + y*y + z*z);
cross(h[0],0.0,0.0,h[4],h[3],h[2],x,y,z);
area[1] = sqrt(x*x + y*y + z*z);
cross(h[5],h[1],0.0,h[4],h[3],h[2],x,y,z);
area[2] = sqrt(x*x + y*y + z*z);
}
double bestsurf = 2.0 * (xprd*yprd + yprd*zprd + zprd*xprd);
double bestsurf = 2.0 * (area[0]+area[1]+area[2]);
// loop thru all possible factorizations of nprocs
// surf = surface area of a proc sub-domain
// for 2d, insure ipz = 1
int ipx,ipy,ipz,nremain;
double surf;
ipx = 1;
while (ipx <= nprocs) {
if (nprocs % ipx == 0) {
@ -796,10 +867,7 @@ void Comm::procs2box()
if (nremain % ipy == 0) {
ipz = nremain/ipy;
if (force->dimension == 3 || ipz == 1) {
boxx = xprd/ipx;
boxy = yprd/ipy;
boxz = zprd/ipz;
surf = boxx*boxy + boxy*boxz + boxz*boxx;
surf = area[0]/ipx/ipy + area[1]/ipx/ipz + area[2]/ipy/ipz;
if (surf < bestsurf) {
bestsurf = surf;
procgrid[0] = ipx;
@ -815,6 +883,19 @@ void Comm::procs2box()
}
}
/* ----------------------------------------------------------------------
vector cross product: c = a x b
------------------------------------------------------------------------- */
void Comm::cross(double ax, double ay, double az,
double bx, double by, double bz,
double &cx, double &cy, double &cz)
{
cx = ay*bz - az*by;
cy = az*bx - ax*bz;
cz = ax*by - ay*bx;
}
/* ----------------------------------------------------------------------
realloc the size of the send buffer as needed with BUFFACTOR & BUFEXTRA
if flag = 1, realloc
@ -896,7 +977,10 @@ void Comm::allocate_swap(int n)
slablo = (double *) memory->smalloc(n*sizeof(double),"comm:slablo");
slabhi = (double *) memory->smalloc(n*sizeof(double),"comm:slabhi");
firstrecv = (int *) memory->smalloc(n*sizeof(int),"comm:firstrecv");
pbc_flags = (int **) memory->create_2d_int_array(n,4,"comm:pbc_flags");
pbc_flag = (int *) memory->smalloc(n*sizeof(int),"comm:pbc_flag");
pbc_dist = (double **) memory->create_2d_double_array(n,3,"comm:pbc_dist");
pbc_dist_border = (double **)
memory->create_2d_double_array(n,3,"comm:pbc_dist_border");
}
/* ----------------------------------------------------------------------
@ -915,7 +999,9 @@ void Comm::free_swap()
memory->sfree(slablo);
memory->sfree(slabhi);
memory->sfree(firstrecv);
memory->destroy_2d_int_array(pbc_flags);
memory->sfree(pbc_flag);
memory->destroy_2d_double_array(pbc_dist);
memory->destroy_2d_double_array(pbc_dist_border);
}
/* ----------------------------------------------------------------------

10
src/comm.h
View File

@ -51,6 +51,7 @@ class Comm : protected Pointers {
int memory_usage(); // tally memory usage
private:
int triclinic; // 0 if domain is orthog, 1 if triclinic
int maxswap; // max # of swaps memory is allocated for
int *sendnum,*recvnum; // # of atoms to send/recv in each swap
int *sendproc,*recvproc; // proc to send/recv to/from at each swap
@ -58,9 +59,9 @@ class Comm : protected Pointers {
int *size_reverse_send; // # to send in each reverse comm
int *size_reverse_recv; // # to recv in each reverse comm
double *slablo,*slabhi; // bounds of slab to send at each swap
int **pbc_flags; // flags for sending atoms thru PBC
// [0] = 1 if any dim is across PBC
// [123] = 1 if dim 012 is across PBC
int *pbc_flag; // flags for sending atoms thru PBC
double **pbc_dist; // distance to adjust atoms coords for PBC
double **pbc_dist_border; // PBC distance for borders()
int comm_x_only,comm_f_only; // 1 if only exchange x,f in for/rev comm
int map_style; // non-0 if global->local mapping is done
@ -74,6 +75,9 @@ class Comm : protected Pointers {
int maxforward,maxreverse; // max # of datums in forward/reverse comm
void procs2box(); // map procs to 3d box
void cross(double, double, double,
double, double, double,
double &, double &, double &); // cross product
void grow_send(int,int); // reallocate send buffer
void grow_recv(int); // free/allocate recv buffer
void grow_list(int, int); // reallocate one sendlist

2
src/compute.h
View File

@ -55,7 +55,7 @@ class Compute : protected Pointers {
virtual void compute_vector() {}
virtual void compute_peratom() {}
virtual int pack_comm(int, int *, double *, int *) {return 0;}
virtual int pack_comm(int, int *, double *, int, double *) {return 0;}
virtual void unpack_comm(int, int, double *) {}
virtual int pack_reverse_comm(int, int, double *) {return 0;}
virtual void unpack_reverse_comm(int, int *, double *) {}

158
src/create_atoms.cpp
View File

@ -54,8 +54,8 @@ void CreateAtoms::command(int narg, char **arg)
error->all("Invalid atom type in create_atoms command");
for (int i = 0; i < nbasis; i++) basistype[i] = itype;
regionflag = -1;
nhybrid = 0;
int regionflag = -1;
int nhybrid = 0;
int iarg = 1;
while (iarg < narg) {
@ -93,30 +93,42 @@ void CreateAtoms::command(int narg, char **arg)
// convert 8 corners of my sub-box from box coords to lattice coords
// min to max = bounding box around the pts in lattice space
subxlo = domain->subxlo;
subxhi = domain->subxhi;
subylo = domain->subylo;
subyhi = domain->subyhi;
subzlo = domain->subzlo;
subzhi = domain->subzhi;
int triclinic = domain->triclinic;
double *bboxlo,*bboxhi;
if (triclinic == 0) {
bboxlo = domain->sublo;
bboxhi = domain->subhi;
} else {
bboxlo = domain->sublo_bound;
bboxhi = domain->subhi_bound;
}
double xmin,ymin,zmin,xmax,ymax,zmax;
xmin = ymin = zmin = BIG;
xmax = ymax = zmax = -BIG;
domain->lattice->bbox(1,subxlo,subylo,subzlo,xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,subxhi,subylo,subzlo,xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,subxlo,subyhi,subzlo,xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,subxhi,subyhi,subzlo,xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,subxlo,subylo,subzhi,xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,subxhi,subylo,subzhi,xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,subxlo,subyhi,subzhi,xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,subxhi,subyhi,subzhi,xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxlo[2],
xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxlo[2],
xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxlo[2],
xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxlo[2],
xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxhi[2],
xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxhi[2],
xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxhi[2],
xmin,ymin,zmin,xmax,ymax,zmax);
domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxhi[2],
xmin,ymin,zmin,xmax,ymax,zmax);
// ilo:ihi,jlo:jhi,klo:khi = loop bounds for lattice overlap of my sub-box
// overlap = any part of a unit cell (face,edge,pt) in common with my sub-box
// in lattice space, sub-box is a tilted box
// but bbox of sub-box is aligned with lattice axes
// ilo:ihi,jlo:jhi,klo:khi = loop bounds for lattice overlap of my subbox
// overlap = any part of a unit cell (face,edge,pt) in common with my subbox
// in lattice space, subbox is a tilted box
// but bbox of subbox is aligned with lattice axes
// so ilo:khi unit cells should completely tile bounding box
// decrement lo values if min < 0, since static_cast(-1.5) = -1
@ -132,25 +144,84 @@ void CreateAtoms::command(int narg, char **arg)
if (ymin < 0.0) jlo--;
if (zmin < 0.0) klo--;
// set bounds for my proc
// if periodic and I am lo/hi proc, adjust bounds by EPSILON
// on lower boundary, allows triclinic atoms just outside box to be added
// on upper boundary, prevents atoms with lower images from being added
double sublo[3],subhi[3];
if (triclinic == 0) {
sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0];
sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1];
sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2];
} else {
sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0];
sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1];
sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2];
}
if (domain->xperiodic) {
if (triclinic && comm->myloc[0] == 0) sublo[0] -= EPSILON;
if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] -= EPSILON;
}
if (domain->yperiodic) {
if (triclinic && comm->myloc[1] == 0) sublo[1] -= EPSILON;
if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] -= EPSILON;
}
if (domain->zperiodic) {
if (triclinic && comm->myloc[2] == 0) sublo[2] -= EPSILON;
if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] -= EPSILON;
}
// iterate on 3d periodic lattice using loop bounds
// invoke add_atom for nbasis atoms in each unit cell
// add_atom converts lattice coords to box coords, checks if in my sub-box
boxxhi = domain->boxxhi;
boxyhi = domain->boxyhi;
boxzhi = domain->boxzhi;
// converts lattice coords to box coords
// add atom if it meets all criteria
double natoms_previous = atom->natoms;
int nlocal_previous = atom->nlocal;
double **basis = domain->lattice->basis;
double x[3],lamda[3];
double *coord;
int i,j,k,m;
for (k = klo; k <= khi; k++)
for (j = jlo; j <= jhi; j++)
for (i = ilo; i <= ihi; i++)
for (m = 0; m < nbasis; m++)
add_atom(basistype[m],i+basis[m][0],j+basis[m][1],k+basis[m][2]);
for (m = 0; m < nbasis; m++) {
x[0] = i + basis[m][0];
x[1] = j + basis[m][1];
x[2] = k + basis[m][2];
// convert from lattice coords to box coords
domain->lattice->lattice2box(x[0],x[1],x[2]);
// if a region was specified, test if atom is in it
if (regionflag >= 0)
if (!domain->regions[regionflag]->match(x[0],x[1],x[2])) continue;
// test if atom is in my subbox
if (triclinic) {
domain->x2lamda(x,lamda);
coord = lamda;
} else coord = x;
if (coord[0] < sublo[0] || coord[0] >= subhi[0] ||
coord[1] < sublo[1] || coord[1] >= subhi[1] ||
coord[2] < sublo[2] || coord[2] >= subhi[2]) continue;
// add the atom to my list of atoms
atom->avec->create_atom(basistype[m],x,nhybrid);
}
// clean up
delete [] basistype;
@ -189,36 +260,3 @@ void CreateAtoms::command(int narg, char **arg)
}
}
}
/* ----------------------------------------------------------------------
add an atom of type at lattice coords x,y,z if it meets all criteria
------------------------------------------------------------------------- */
void CreateAtoms::add_atom(int ntype, double x, double y, double z)
{
// convert from lattice coords to box coords
domain->lattice->lattice2box(x,y,z);
// if a region was specified, test if atom is in it
if (regionflag >= 0)
if (!domain->regions[regionflag]->match(x,y,z)) return;
// test if atom is in my subbox
if (x < subxlo || x >= subxhi ||
y < subylo || y >= subyhi ||
z < subzlo || z >= subzhi) return;
// don't put atoms within EPSILON of upper periodic boundary
// else may overlap image atom at lower boundary
if (domain->xperiodic && fabs(x-boxxhi) < EPSILON) return;
if (domain->yperiodic && fabs(y-boxyhi) < EPSILON) return;
if (domain->zperiodic && fabs(z-boxzhi) < EPSILON) return;
// add the atom to my list of atoms
atom->avec->create_atom(ntype,x,y,z,nhybrid);
}

7
src/create_atoms.h
View File

@ -22,13 +22,6 @@ class CreateAtoms : protected Pointers {
public:
CreateAtoms(class LAMMPS *);
void command(int, char **);
private:
int regionflag,nhybrid;
double boxxhi,boxyhi,boxzhi;
double subxlo,subxhi,subylo,subyhi,subzlo,subzhi;
void add_atom(int, double, double, double);
};
}

47
src/create_box.cpp
View File

@ -18,6 +18,7 @@
#include "force.h"
#include "domain.h"
#include "region.h"
#include "region_prism.h"
#include "comm.h"
#include "update.h"
#include "error.h"
@ -39,6 +40,8 @@ void CreateBox::command(int narg, char **arg)
if (force->dimension == 2 && domain->zperiodic == 0)
error->all("Cannot run 2d simulation with nonperiodic Z dimension");
domain->box_exist = 1;
// find region ID
int iregion;
@ -50,26 +53,33 @@ void CreateBox::command(int narg, char **arg)
if (domain->regions[iregion]->interior == 0)
error->all("Create_box region must be of type inside");
// set global box from region extent
// if region not prism:
// setup orthogonal domain
// set simulation domain from region extent
// if region is prism:
// seutp triclinic domain
// set simulation domain params from prism params
domain->boxxlo = domain->regions[iregion]->extent_xlo;
domain->boxxhi = domain->regions[iregion]->extent_xhi;
domain->boxylo = domain->regions[iregion]->extent_ylo;
domain->boxyhi = domain->regions[iregion]->extent_yhi;
domain->boxzlo = domain->regions[iregion]->extent_zlo;
domain->boxzhi = domain->regions[iregion]->extent_zhi;
domain->box_exist = 1;
if (strcmp(domain->regions[iregion]->style,"prism") != 0) {
domain->boxxlo = domain->regions[iregion]->extent_xlo;
domain->boxxhi = domain->regions[iregion]->extent_xhi;
domain->boxylo = domain->regions[iregion]->extent_ylo;
domain->boxyhi = domain->regions[iregion]->extent_yhi;
domain->boxzlo = domain->regions[iregion]->extent_zlo;
domain->boxzhi = domain->regions[iregion]->extent_zhi;
if (comm->me == 0) {
if (screen)
fprintf(screen,"Created box = (%g %g %g) to (%g %g %g)\n",
domain->boxxlo,domain->boxylo,domain->boxzlo,
domain->boxxhi,domain->boxyhi,domain->boxzhi);
if (logfile)
fprintf(logfile,"Created box = (%g %g %g) to (%g %g %g)\n",
domain->boxxlo,domain->boxylo,domain->boxzlo,
domain->boxxhi,domain->boxyhi,domain->boxzhi);
} else {
domain->triclinic = 1;
RegPrism *region = (RegPrism *) domain->regions[iregion];
domain->boxxlo = region->xlo;
domain->boxxhi = region->xhi;
domain->boxylo = region->ylo;
domain->boxyhi = region->yhi;
domain->boxzlo = region->zlo;
domain->boxzhi = region->zhi;
domain->xy = region->xy;
domain->xz = region->xz;
domain->yz = region->yz;
}
// if molecular, zero out topology info
@ -100,6 +110,7 @@ void CreateBox::command(int narg, char **arg)
atom->allocate_type_arrays();
domain->print_box("Created ");
domain->set_initial_box();
domain->set_global_box();
comm->set_procs();

4
src/delete_atoms.cpp
View File

@ -196,13 +196,15 @@ void DeleteAtoms::delete_overlap(int narg, char **arg, int *list)
// setup domain, communication and neighboring
// acquire ghosts
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
// call to build() forces memory allocation for neighbor lists
// build half list explicitly if build() doesn't do it

2
src/delete_bonds.cpp
View File

@ -104,11 +104,13 @@ void DeleteBonds::command(int narg, char **arg)
// border swap to insure type and mask is current for off-proc atoms
// enforce PBC before in case atoms are outside box
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
// set topology interactions either off or on
// criteria for an interaction to potentially be changed (set flag = 1)

2
src/displace_atoms.cpp
View File

@ -165,10 +165,12 @@ void DisplaceAtoms::command(int narg, char **arg)
// move atoms to new processors
// enforce PBC before in case atoms are outside box
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
comm->exchange();
if (domain->triclinic) domain->lamda2x(atom->nlocal);
// check if any atoms were lost

522
src/domain.cpp
View File

@ -58,8 +58,14 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
boundary[1][0] = boundary[1][1] = 0;
boundary[2][0] = boundary[2][1] = 0;
triclinic = 0;
boxxlo = boxylo = boxzlo = -0.5;
boxxhi = boxyhi = boxzhi = 0.5;
xy = xz = yz = 0.0;
prd_lamda[0] = prd_lamda[1] = prd_lamda[2] = 1.0;
boxlo_lamda[0] = boxlo_lamda[1] = boxlo_lamda[2] = 0.0;
boxhi_lamda[0] = boxhi_lamda[1] = boxhi_lamda[2] = 1.0;
lattice = NULL;
nregion = maxregion = 0;
@ -79,26 +85,47 @@ Domain::~Domain()
void Domain::init()
{
// set box_change if box dimensions ever change
// due to shrink-wrapping, fix nph, npt, volume/rescale, uniaxial
// set box_change if box dimensions/shape ever changes
// due to shrink-wrapping, fixes that change volume (npt, vol/rescale, etc)
box_change = 0;
if (nonperiodic == 2) box_change = 1;
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"nph") == 0) box_change = 1;
if (strcmp(modify->fix[i]->style,"npt") == 0) box_change = 1;
if (strcmp(modify->fix[i]->style,"volume/rescale") == 0) box_change = 1;
if (strcmp(modify->fix[i]->style,"uniaxial") == 0) box_change = 1;
}
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->box_change) box_change = 1;
}
/* ----------------------------------------------------------------------
set initial global box from boxlo/hi (already set by caller)
adjust for shrink-wrapped dims
set initial global box
assumes boxlo/hi and triclinic tilts are already set
------------------------------------------------------------------------- */
void Domain::set_initial_box()
{
// error checks for orthogonal and triclinic domains
if (boxxlo >= boxxhi || boxylo >= boxyhi || boxzlo >= boxzhi)
error->one("Box bounds are invalid");
if (triclinic) {
if (force->dimension == 2 && (xz != 0.0 || yz != 0.0))
error->all("Cannot skew triclinic box in z for 2d simulation");
if (xy != 0.0 && (!xperiodic || !yperiodic))
error->all("Triclinic box must be periodic in skewed dimensions");
if (xz != 0.0 && (!xperiodic || !zperiodic))
error->all("Triclinic box must be periodic in skewed dimensions");
if (yz != 0.0 && (!yperiodic || !zperiodic))
error->all("Triclinic box must be periodic in skewed dimensions");
if (fabs(xy/(boxyhi-boxylo)) > 0.5)
error->all("Triclinic box skew is too large");
if (fabs(xz/(boxzhi-boxzlo)) > 0.5)
error->all("Triclinic box skew is too large");
if (fabs(yz/(boxzhi-boxzlo)) > 0.5)
error->all("Triclinic box skew is too large");
}
// adjust box lo/hi for shrink-wrapped dims
if (boundary[0][0] == 2) boxxlo -= SMALL;
else if (boundary[0][0] == 3) minxlo = boxxlo;
if (boundary[0][1] == 2) boxxhi += SMALL;
@ -116,54 +143,132 @@ void Domain::set_initial_box()
}
/* ----------------------------------------------------------------------
set global prd, prd_half, prd[], boxlo/hi[] from boxlo/hi
set global box params
assumes boxlo/hi and triclinic tilts are already set
------------------------------------------------------------------------- */
void Domain::set_global_box()
{
xprd = boxxhi - boxxlo;
yprd = boxyhi - boxylo;
zprd = boxzhi - boxzlo;
prd[0] = xprd = boxxhi - boxxlo;
prd[1] = yprd = boxyhi - boxylo;
prd[2] = zprd = boxzhi - boxzlo;
xprd_half = 0.5*xprd;
yprd_half = 0.5*yprd;
zprd_half = 0.5*zprd;
prd[0] = xprd; prd[1] = yprd; prd[2] = zprd;
boxlo[0] = boxxlo; boxlo[1] = boxylo; boxlo[2] = boxzlo;
boxhi[0] = boxxhi; boxhi[1] = boxyhi; boxhi[2] = boxzhi;
if (triclinic) {
h[0] = xprd;
h[1] = yprd;
h[2] = zprd;
h[3] = yz;
h[4] = xz;
h[5] = xy;
h_inv[0] = 1.0/h[0];
h_inv[1] = 1.0/h[1];
h_inv[2] = 1.0/h[2];
h_inv[3] = -h[3] / (h[1]*h[2]);
h_inv[4] = (h[3]*h[5] - h[1]*h[4]) / (h[0]*h[1]*h[2]);
h_inv[5] = -h[5] / (h[0]*h[1]);
boxlo_bound[0] = MIN(boxlo[0],boxlo[0]+xy);
boxlo_bound[0] = MIN(boxlo_bound[0],boxlo_bound[0]+xz);
boxlo_bound[1] = MIN(boxlo[1],boxlo[1]+yz);
boxlo_bound[2] = boxlo[2];
boxhi_bound[0] = MAX(boxhi[0],boxhi[0]+xy);
boxhi_bound[0] = MAX(boxhi_bound[0],boxhi_bound[0]+xz);
boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz);
boxhi_bound[2] = boxhi[2];
}
}
/* ----------------------------------------------------------------------
set local subxlo-subzhi, sublo/hi[] from global boxxlo-boxzhi and proc grid
set lamda box params, only need be done one time
assumes global box is defined and proc assignment has been made by comm
for uppermost proc, insure subhi = 1.0 (in case round-off occurs)
------------------------------------------------------------------------- */
void Domain::set_lamda_box()
{
int *myloc = comm->myloc;
int *procgrid = comm->procgrid;
sublo_lamda[0] = 1.0*myloc[0] / procgrid[0];
sublo_lamda[1] = 1.0*myloc[1] / procgrid[1];
sublo_lamda[2] = 1.0*myloc[2] / procgrid[2];
if (myloc[0] < procgrid[0]-1)
subhi_lamda[0] = 1.0*(myloc[0]+1) / procgrid[0];
else subhi_lamda[0] = 1.0;
if (myloc[1] < procgrid[1]-1)
subhi_lamda[1] = 1.0*(myloc[1]+1) / procgrid[1];
else subhi_lamda[1] = 1.0;
if (myloc[2] < procgrid[2]-1)
subhi_lamda[2] = 1.0*(myloc[2]+1) / procgrid[2];
else subhi_lamda[2] = 1.0;
}
/* ----------------------------------------------------------------------
set local subbox params
assumes global box is defined and proc assignment has been made
for uppermost proc, insure subhi = boxhi (in case round-off occurs)
for triclinic, lo/hi_bound is a bbox around all 8 tilted pts of sub-box
------------------------------------------------------------------------- */
void Domain::set_local_box()
{
subxlo = boxxlo + comm->myloc[0] * xprd / comm->procgrid[0];
if (comm->myloc[0] < comm->procgrid[0]-1)
subxhi = boxxlo + (comm->myloc[0]+1) * xprd / comm->procgrid[0];
else subxhi = boxxhi;
int *myloc = comm->myloc;
int *procgrid = comm->procgrid;
subylo = boxylo + comm->myloc[1] * yprd / comm->procgrid[1];
if (comm->myloc[1] < comm->procgrid[1]-1)
subyhi = boxylo + (comm->myloc[1]+1) * yprd / comm->procgrid[1];
else subyhi = boxyhi;
if (domain->triclinic == 0) {
sublo[0] = boxlo[0] + myloc[0] * xprd / procgrid[0];
if (myloc[0] < procgrid[0]-1)
subhi[0] = boxlo[0] + (myloc[0]+1) * xprd / procgrid[0];
else subhi[0] = boxhi[0];
sublo[1] = boxlo[1] + myloc[1] * yprd / procgrid[1];
if (myloc[1] < procgrid[1]-1)
subhi[1] = boxlo[1] + (myloc[1]+1) * yprd / procgrid[1];
else subhi[1] = boxhi[1];
subzlo = boxzlo + comm->myloc[2] * zprd / comm->procgrid[2];
if (comm->myloc[2] < comm->procgrid[2]-1)
subzhi = boxzlo + (comm->myloc[2]+1) * zprd / comm->procgrid[2];
else subzhi = boxzhi;
sublo[2] = boxlo[2] + myloc[2] * zprd / procgrid[2];
if (myloc[2] < procgrid[2]-1)
subhi[2] = boxlo[2] + (myloc[2]+1) * zprd / procgrid[2];
else subhi[2] = boxhi[2];
sublo[0] = subxlo; sublo[1] = subylo; sublo[2] = subzlo;
subhi[0] = subxhi; subhi[1] = subyhi; subhi[2] = subzhi;
} else {
sublo[0] = h[0]*sublo_lamda[0] + h[5]*sublo_lamda[1] +
h[4]*sublo_lamda[2] + boxlo[0];
sublo[1] = h[1]*sublo_lamda[1] + h[3]*sublo_lamda[2] + boxlo[1];
sublo[2] = h[2]*sublo_lamda[2] + boxlo[2];
subhi[0] = sublo[0] + xprd/procgrid[0];
subhi[1] = sublo[1] + yprd/procgrid[1];
subhi[2] = sublo[2] + zprd/procgrid[2];
sublo_bound[0] = MIN(sublo[0],sublo[0] + xy/procgrid[1]);
sublo_bound[0] = MIN(sublo_bound[0],sublo_bound[0] + xz/procgrid[2]);
sublo_bound[1] = MIN(sublo[1],sublo[1] + yz/procgrid[2]);
sublo_bound[2] = sublo[2];
subhi_bound[0] = MAX(subhi[0],subhi[0] + xy/procgrid[1]);
subhi_bound[0] = MAX(subhi_bound[0],subhi_bound[0] + xz/procgrid[2]);
subhi_bound[1] = MAX(subhi[1],subhi[1] + yz/procgrid[2]);
subhi_bound[2] = subhi[2];
}
}
/* ----------------------------------------------------------------------
reset global & local boxes due to global box boundary changes
if shrink-wrapped, determine atom extent and reset boxlo/hi
set global & local boxes from new boxlo/hi values
shrink-wrapping only occurs in non-periodic, non-triclinic dims
------------------------------------------------------------------------- */
void Domain::reset_box()
@ -199,6 +304,7 @@ void Domain::reset_box()
MPI_Allreduce(extent,all,6,MPI_DOUBLE,MPI_MAX,world);
// in shrink-wrapped dims, set box by atom extent
// if set, observe min box size settings
if (xperiodic == 0) {
if (boundary[0][0] == 2) boxxlo = -all[0][0] - SMALL;
@ -229,9 +335,10 @@ void Domain::reset_box()
/* ----------------------------------------------------------------------
enforce PBC and modify box image flags for each atom
called every reneighboring
called every reneighboring and by other commands that change atoms
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
for triclinic, atoms must be in lamda coords (0-1) before pbc is called
image = 10 bits for each dimension
increment/decrement in wrap-around fashion
------------------------------------------------------------------------- */
@ -239,23 +346,34 @@ void Domain::reset_box()
void Domain::pbc()
{
int i,idim,otherdims;
double *lo,*hi,*period;
int nlocal = atom->nlocal;
double **x = atom->x;
int *image = atom->image;
if (xperiodic) {
for (i = 0; i < nlocal; i++) {
if (x[i][0] < boxxlo) {
x[i][0] += xprd;
if (triclinic == 0) {
lo = boxlo;
hi = boxhi;
period = prd;
} else {
lo = boxlo_lamda;
hi = boxhi_lamda;
period = prd_lamda;
}
for (i = 0; i < nlocal; i++) {
if (xperiodic) {
if (x[i][0] < lo[0]) {
x[i][0] += period[0];
idim = image[i] & 1023;
otherdims = image[i] ^ idim;
idim--;
idim &= 1023;
image[i] = otherdims | idim;
}
if (x[i][0] >= boxxhi) {
x[i][0] -= xprd;
x[i][0] = MAX(x[i][0],boxxlo);
if (x[i][0] >= hi[0]) {
x[i][0] -= period[0];
x[i][0] = MAX(x[i][0],lo[0]);
idim = image[i] & 1023;
otherdims = image[i] ^ idim;
idim++;
@ -263,21 +381,19 @@ void Domain::pbc()
image[i] = otherdims | idim;
}
}
}
if (yperiodic) {
for (i = 0; i < nlocal; i++) {
if (x[i][1] < boxylo) {
x[i][1] += yprd;
if (yperiodic) {
if (x[i][1] < lo[1]) {
x[i][1] += period[1];
idim = (image[i] >> 10) & 1023;
otherdims = image[i] ^ (idim << 10);
idim--;
idim &= 1023;
image[i] = otherdims | (idim << 10);
}
if (x[i][1] >= boxyhi) {
x[i][1] -= yprd;
x[i][1] = MAX(x[i][1],boxylo);
if (x[i][1] >= hi[1]) {
x[i][1] -= period[1];
x[i][1] = MAX(x[i][1],lo[1]);
idim = (image[i] >> 10) & 1023;
otherdims = image[i] ^ (idim << 10);
idim++;
@ -285,21 +401,19 @@ void Domain::pbc()
image[i] = otherdims | (idim << 10);
}
}
}
if (zperiodic) {
for (i = 0; i < nlocal; i++) {
if (x[i][2] < boxzlo) {
x[i][2] += zprd;
if (zperiodic) {
if (x[i][2] < lo[2]) {
x[i][2] += period[2];
idim = image[i] >> 20;
otherdims = image[i] ^ (idim << 20);
idim--;
idim &= 1023;
image[i] = otherdims | (idim << 20);
}
if (x[i][2] >= boxzhi) {
x[i][2] -= zprd;
x[i][2] = MAX(x[i][2],boxzlo);
if (x[i][2] >= hi[2]) {
x[i][2] -= period[2];
x[i][2] = MAX(x[i][2],lo[2]);
idim = image[i] >> 20;
otherdims = image[i] ^ (idim << 20);
idim++;
@ -313,53 +427,123 @@ void Domain::pbc()
/* ----------------------------------------------------------------------
minimum image convention
use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed
------------------------------------------------------------------------- */
void Domain::minimum_image(double *dx, double *dy, double *dz)
void Domain::minimum_image(double &dx, double &dy, double &dz)
{
if (xperiodic) {
if (fabs(*dx) > xprd_half) {
if (*dx < 0.0) *dx += xprd;
else *dx -= xprd;
if (triclinic == 0) {
if (xperiodic) {
if (fabs(dx) > xprd_half) {
if (dx < 0.0) dx += xprd;
else dx -= xprd;
}
}
}
if (yperiodic) {
if (fabs(*dy) > yprd_half) {
if (*dy < 0.0) *dy += yprd;
else *dy -= yprd;
if (yperiodic) {
if (fabs(dy) > yprd_half) {
if (dy < 0.0) dy += yprd;
else dy -= yprd;
}
}
}
if (zperiodic) {
if (fabs(*dz) > zprd_half) {
if (*dz < 0.0) *dz += zprd;
else *dz -= zprd;
if (zperiodic) {
if (fabs(dz) > zprd_half) {
if (dz < 0.0) dz += zprd;
else dz -= zprd;
}
}
} else {
if (zperiodic) {
if (fabs(dz) > zprd_half) {
if (dz < 0.0) {
dz += zprd;
dy += yz;
dx += xz;
} else {
dz -= zprd;
dy -= yz;
dx -= xz;
}
}
}
if (yperiodic) {
if (fabs(dy) > yprd_half) {
if (dy < 0.0) {
dy += yprd;
dx += xy;
} else {
dy -= yprd;
dx -= xy;
}
}
}
if (xperiodic) {
if (fabs(dx) > xprd_half) {
if (dx < 0.0) dx += xprd;
else dx -= xprd;
}
}
}
}
/* ----------------------------------------------------------------------
minimum image convention
use 1/2 of box size as test
use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed
------------------------------------------------------------------------- */
void Domain::minimum_image(double *x)
void Domain::minimum_image(double *delta)
{
if (xperiodic) {
if (fabs(x[0]) > xprd_half) {
if (x[0] < 0.0) x[0] += xprd;
else x[0] -= xprd;
if (triclinic == 0) {
if (xperiodic) {
if (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
}
if (yperiodic) {
if (fabs(x[1]) > yprd_half) {
if (x[1] < 0.0) x[1] += yprd;
else x[1] -= yprd;
if (yperiodic) {
if (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) delta[1] += yprd;
else delta[1] -= yprd;
}
}
}
if (zperiodic) {
if (fabs(x[2]) > zprd_half) {
if (x[2] < 0.0) x[2] += zprd;
else x[2] -= zprd;
if (zperiodic) {
if (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) delta[2] += zprd;
else delta[2] -= zprd;
}
}
} else {
if (zperiodic) {
if (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) {
delta[2] += zprd;
delta[1] += yz;
delta[0] += xz;
} else {
delta[2] -= zprd;
delta[1] -= yz;
delta[0] -= xz;
}
}
}
if (yperiodic) {
if (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) {
delta[1] += yprd;
delta[0] += xy;
} else {
delta[1] -= yprd;
delta[0] -= xy;
}
}
}
if (xperiodic) {
if (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
}
}
@ -369,87 +553,113 @@ void Domain::minimum_image(double *x)
adjust image accordingly
resulting coord must satisfy lo <= coord < hi
MAX is important since coord - prd < lo can happen when coord = hi
for triclinic, convert atom to lamda coords (0-1) before doing remap
image = 10 bits for each dimension
increment/decrement in wrap-around fashion
------------------------------------------------------------------------- */
void Domain::remap(double &x, double &y, double &z, int &image)
void Domain::remap(double *x, int &image)
{
double *lo,*hi,*period,*coord;
double lamda[3];
if (triclinic == 0) {
lo = boxlo;
hi = boxhi;
period = prd;
coord = x;
} else {
lo = boxlo_lamda;
hi = boxhi_lamda;
period = prd_lamda;
x2lamda(x,lamda);
coord = lamda;
}
if (xperiodic) {
while (x < boxxlo) {
x += xprd;
while (coord[0] < lo[0]) {
coord[0] += period[0];
int idim = image & 1023;
int otherdims = image ^ idim;
idim--;
idim &= 1023;
image = otherdims | idim;
}
while (x >= boxxhi) {
x -= xprd;
while (coord[0] >= hi[0]) {
coord[0] -= period[0];
int idim = image & 1023;
int otherdims = image ^ idim;
idim++;
idim &= 1023;
image = otherdims | idim;
}
x = MAX(x,boxxlo);
coord[0] = MAX(coord[0],lo[0]);
}
if (yperiodic) {
while (y < boxylo) {
y += yprd;
while (coord[1] < lo[1]) {
coord[1] += period[1];
int idim = (image >> 10) & 1023;
int otherdims = image ^ (idim << 10);
idim--;
idim &= 1023;
image = otherdims | (idim << 10);
}
while (y >= boxyhi) {
y -= yprd;
while (coord[1] >= hi[1]) {
coord[1] -= period[1];
int idim = (image >> 10) & 1023;
int otherdims = image ^ (idim << 10);
idim++;
idim &= 1023;
image = otherdims | (idim << 10);
}
y = MAX(y,boxylo);
coord[1] = MAX(coord[1],lo[1]);
}
if (zperiodic) {
while (z < boxzlo) {
z += zprd;
while (coord[2] < lo[2]) {
coord[2] += period[2];
int idim = image >> 20;
int otherdims = image ^ (idim << 20);
idim--;
idim &= 1023;
image = otherdims | (idim << 20);
}
while (z >= boxzhi) {
z -= zprd;
while (coord[2] >= hi[2]) {
coord[2] -= period[2];
int idim = image >> 20;
int otherdims = image ^ (idim << 20);
idim++;
idim &= 1023;
image = otherdims | (idim << 20);
}
z = MAX(z,boxzlo);
coord[2] = MAX(coord[2],lo[2]);
}
if (triclinic) lamda2x(coord,x);
}
/* ----------------------------------------------------------------------
unmap the point via image flags
don't reset image flag
for triclinic, use h[] to add in tilt factors in other dims as needed
------------------------------------------------------------------------- */
void Domain::unmap(double &x, double &y, double &z, int image)
void Domain::unmap(double *x, int image)
{
int xbox = (image & 1023) - 512;
int ybox = (image >> 10 & 1023) - 512;
int zbox = (image >> 20) - 512;
x = x + xbox*xprd;
y = y + ybox*yprd;
z = z + zbox*zprd;
if (triclinic == 0) {
x[0] += xbox*xprd;
x[1] += ybox*yprd;
x[2] += zbox*zprd;
} else {
x[0] += h[0]*xbox + h[5]*ybox + h[4]*zbox;
x[1] += h[1]*ybox + h[3]*zbox;
x[2] += h[2]*zbox;
}
}
/* ----------------------------------------------------------------------
@ -551,3 +761,103 @@ void Domain::set_boundary(int narg, char **arg)
boundary[2][0] >= 2 || boundary[2][1] >= 2) nonperiodic = 2;
}
}
/* ----------------------------------------------------------------------
print box info, orthogonal or triclinic
------------------------------------------------------------------------- */
void Domain::print_box(char *str)
{
if (comm->me == 0) {
if (screen) {
if (domain->triclinic == 0)
fprintf(screen,"%sorthogonal box = (%g %g %g) to (%g %g %g)\n",
str,boxxlo,boxylo,boxzlo,boxxhi,boxyhi,boxzhi);
else {
char *format =
"%striclinic box = (%g %g %g) to (%g %g %g) with tilt (%g %g %g)\n";
fprintf(screen,format,
str,boxxlo,boxylo,boxzlo,boxxhi,boxyhi,boxzhi,xy,xz,yz);
}
}
if (logfile) {
if (triclinic == 0)
fprintf(logfile,"%sorthogonal box = (%g %g %g) to (%g %g %g)\n",
str,boxxlo,boxylo,boxzlo,boxxhi,boxyhi,boxzhi);
else {
char *format =
"%striclinic box = (%g %g %g) to (%g %g %g) with tilt (%g %g %g)\n";
fprintf(logfile,format,
str,boxxlo,boxylo,boxzlo,boxxhi,boxyhi,boxzhi,xy,xz,yz);
}
}
}
}
/* ----------------------------------------------------------------------
convert triclinic 0-1 lamda coords to box coords for all N atoms
x = H lamda + x0;
------------------------------------------------------------------------- */
void Domain::lamda2x(int n)
{
double **x = atom->x;
for (int i = 0; i < n; i++) {
x[i][0] = h[0]*x[i][0] + h[5]*x[i][1] + h[4]*x[i][2] + boxlo[0];
x[i][1] = h[1]*x[i][1] + h[3]*x[i][2] + boxlo[1];
x[i][2] = h[2]*x[i][2] + boxlo[2];
}
}
/* ----------------------------------------------------------------------
convert box coords to triclinic 0-1 lamda coords for all N atoms
lamda = H^-1 (x - x0)
------------------------------------------------------------------------- */
void Domain::x2lamda(int n)
{
double delta[3];
double **x = atom->x;
for (int i = 0; i < n; i++) {
delta[0] = x[i][0] - boxlo[0];
delta[1] = x[i][1] - boxlo[1];
delta[2] = x[i][2] - boxlo[2];
x[i][0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2];
x[i][1] = h_inv[1]*delta[1] + h_inv[3]*delta[2];
x[i][2] = h_inv[2]*delta[2];
}
}
/* ----------------------------------------------------------------------
convert triclinic 0-1 lamda coords to box coords for one atom
x = H lamda + x0;
lamda and x can point to same 3-vector
------------------------------------------------------------------------- */
void Domain::lamda2x(double *lamda, double *x)
{
x[0] = h[0]*lamda[0] + h[5]*lamda[1] + h[4]*lamda[2] + boxlo[0];
x[1] = h[1]*lamda[1] + h[3]*lamda[2] + boxlo[1];
x[2] = h[2]*lamda[2] + boxlo[2];
}
/* ----------------------------------------------------------------------
convert box coords to triclinic 0-1 lamda coords for one atom
lamda = H^-1 (x - x0)
x and lamda can point to same 3-vector
------------------------------------------------------------------------- */
void Domain::x2lamda(double *x, double *lamda)
{
double delta[3];
delta[0] = x[0] - boxlo[0];
delta[1] = x[1] - boxlo[1];
delta[2] = x[2] - boxlo[2];
lamda[0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2];
lamda[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2];
lamda[2] = h_inv[2]*delta[2];
}

94
src/domain.h
View File

@ -20,61 +20,91 @@ namespace LAMMPS_NS {
class Domain : protected Pointers {
public:
int box_exist; // 0 = not yet created, 1 = exists
int box_exist; // 0 = not yet created, 1 = exists
int nonperiodic; // 0 = periodic in all 3 dims
// 1 = periodic or fixed in all 6
// 2 = shrink-wrap in any of 6
int xperiodic,yperiodic,zperiodic; // 0 = not periodic, 1 = periodic
int boundary[3][2]; // settings for 6 boundaries
// 0 = periodic
// 1 = fixed non-periodic
// 2 = shrink-wrap non-periodic
// 3 = shrink-wrap non-per w/ min
int nonperiodic; // 0 = periodic in all 3 dims
// 1 = periodic or fixed in all 6
// 2 = shrink-wrap in any of 6
int xperiodic,yperiodic,zperiodic; // 0 = non-periodic, 1 = periodic
int periodicity[3]; // xyz periodicity as array
double minxlo,minxhi; // minimum size of global box
double minylo,minyhi; // when shrink-wrapping
double minzlo,minzhi;
int boundary[3][2]; // settings for 6 boundaries
// 0 = periodic
// 1 = fixed non-periodic
// 2 = shrink-wrap non-periodic
// 3 = shrink-wrap non-per w/ min
double boxxlo,boxxhi; // global box boundaries
double boxylo,boxyhi;
int triclinic; // 0 = orthog box, 1 = triclinic
// orthogonal box
double xprd,yprd,zprd; // global box dimensions
double xprd_half,yprd_half,zprd_half; // half dimensions
double prd[3]; // array form of dimensions
// triclinic box
// xprd,half,prd = same
// (as if untilted)
double prd_lamda[3]; // lamda box = (1,1,1)
double boxxlo,boxxhi; // orthogonal box
double boxylo,boxyhi; // global box bounds
double boxzlo,boxzhi;
double boxlo[3],boxhi[3]; // array form of box bounds
double xprd,yprd,zprd; // global box size
double xprd_half,yprd_half,zprd_half;
// triclinic box
// boxxlo/hi,boxlo/hi = same
// (as if untilted)
double boxlo_lamda[3],boxhi_lamda[3]; // lamda box = (0,1)
double boxlo_bound[3],boxhi_bound[3]; // bounding box of tilted domain
double subxlo,subxhi; // sub-box boudaries on this proc
double subylo,subyhi;
double subzlo,subzhi;
// orthogonal box
double minxlo,minxhi; // minimum size of global box
double minylo,minyhi; // when shrink-wrapping
double minzlo,minzhi; // no shrink-wrapping for triclinic
double boxlo[3],boxhi[3]; // global box bounds as arrays
double prd[3]; // global box size as array
double sublo[3],subhi[3]; // sub-box bounds as arrays
int periodicity[3]; // xyz periodic as array
// orthogonal box
double sublo[3],subhi[3]; // sub-box bounds on this proc
int box_change; // 1 if box bounds ever change, 0 if fixed
// triclinic box
// sublo = lower-left corner with tilt
// subhi = upper-right as if untilted
double sublo_lamda[3],subhi_lamda[3]; // bounds of subbox in lamda
double sublo_bound[3],subhi_bound[3]; // bounding box of tilted subdomain
class Lattice *lattice; // user-defined lattice
// triclinic box
double xy,xz,yz; // 3 tilt factors
double h[6],h_inv[6]; // shape matrix in Voigt notation
int nregion; // # of defined Regions
int maxregion; // max # list can hold
class Region **regions; // list of defined Regions
int box_change; // 1 if box bounds ever change, 0 if fixed
class Lattice *lattice; // user-defined lattice
int nregion; // # of defined Regions
int maxregion; // max # list can hold
class Region **regions; // list of defined Regions
Domain(class LAMMPS *);
~Domain();
void init();
void set_initial_box();
void set_global_box();
void set_lamda_box();
void set_local_box();
void reset_box();
void pbc();
void remap(double &, double &, double &, int &);
void unmap(double &, double &, double &, int);
void minimum_image(double *, double *, double *);
void remap(double *, int &);
void unmap(double *, int);
void minimum_image(double &, double &, double &);
void minimum_image(double *);
void set_lattice(int, char **);
void add_region(int, char **);
void set_boundary(int, char **);
void print_box(char *);
void lamda2x(int);
void x2lamda(int);
void lamda2x(double *, double *);
void x2lamda(double *, double *);
};
}

19
src/dump.cpp
View File

@ -18,6 +18,7 @@
#include "dump.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "group.h"
#include "output.h"
#include "memory.h"
@ -118,6 +119,24 @@ void Dump::write()
if (multifile) openfile();
// simulation box bounds
if (domain->triclinic == 0) {
boxxlo = domain->boxxlo;
boxxhi = domain->boxxhi;
boxylo = domain->boxylo;
boxyhi = domain->boxyhi;
boxzlo = domain->boxzlo;
boxzhi = domain->boxzhi;
} else {
boxxlo = domain->boxlo_bound[0];
boxxhi = domain->boxhi_bound[0];
boxylo = domain->boxlo_bound[1];
boxyhi = domain->boxhi_bound[1];
boxzlo = domain->boxlo_bound[2];
boxzhi = domain->boxhi_bound[2];
}
// nme = # of dump lines this proc will contribute to dump
// ntotal = total # of dump lines
// nmax = max # of dump lines on any proc

4
src/dump.h
View File

@ -51,6 +51,10 @@ class Dump : protected Pointers {
virtual int memory_usage();
protected:
double boxxlo,boxxhi; // local copies of domain values
double boxylo,boxyhi; // adjusted to be bounding box for triclinic
double boxzlo,boxzhi;
virtual void openfile();
virtual int modify_param(int, char **) {return 0;}

27
src/dump_atom.cpp
View File

@ -40,6 +40,9 @@ DumpAtom::DumpAtom(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
void DumpAtom::init()
{
if (scale_flag && domain->triclinic)
error->all("Cannot dump scaled coords with triclinic box");
if (image_flag == 0) size_one = 5;
else size_one = 8;
@ -142,12 +145,12 @@ void DumpAtom::header_binary(int ndump)
{
fwrite(&update->ntimestep,sizeof(int),1,fp);
fwrite(&ndump,sizeof(int),1,fp);
fwrite(&domain->boxxlo,sizeof(double),1,fp);
fwrite(&domain->boxxhi,sizeof(double),1,fp);
fwrite(&domain->boxylo,sizeof(double),1,fp);
fwrite(&domain->boxyhi,sizeof(double),1,fp);
fwrite(&domain->boxzlo,sizeof(double),1,fp);
fwrite(&domain->boxzhi,sizeof(double),1,fp);
fwrite(&boxxlo,sizeof(double),1,fp);
fwrite(&boxxhi,sizeof(double),1,fp);
fwrite(&boxylo,sizeof(double),1,fp);
fwrite(&boxyhi,sizeof(double),1,fp);
fwrite(&boxzlo,sizeof(double),1,fp);
fwrite(&boxzhi,sizeof(double),1,fp);
fwrite(&size_one,sizeof(int),1,fp);
if (multiproc) {
int one = 1;
@ -164,9 +167,9 @@ void DumpAtom::header_item(int ndump)
fprintf(fp,"ITEM: NUMBER OF ATOMS\n");
fprintf(fp,"%d\n",ndump);
fprintf(fp,"ITEM: BOX BOUNDS\n");
fprintf(fp,"%g %g\n",domain->boxxlo,domain->boxxhi);
fprintf(fp,"%g %g\n",domain->boxylo,domain->boxyhi);
fprintf(fp,"%g %g\n",domain->boxzlo,domain->boxzhi);
fprintf(fp,"%g %g\n",boxxlo,boxxhi);
fprintf(fp,"%g %g\n",boxylo,boxyhi);
fprintf(fp,"%g %g\n",boxzlo,boxzhi);
fprintf(fp,"ITEM: ATOMS\n");
}
@ -181,9 +184,6 @@ int DumpAtom::pack_scale_image()
double **x = atom->x;
int nlocal = atom->nlocal;
double boxxlo = domain->boxxlo;
double boxylo = domain->boxylo;
double boxzlo = domain->boxzlo;
double invxprd = 1.0/domain->xprd;
double invyprd = 1.0/domain->yprd;
double invzprd = 1.0/domain->zprd;
@ -214,9 +214,6 @@ int DumpAtom::pack_scale_noimage()
double **x = atom->x;
int nlocal = atom->nlocal;
double boxxlo = domain->boxxlo;
double boxylo = domain->boxylo;
double boxzlo = domain->boxzlo;
double invxprd = 1.0/domain->xprd;
double invyprd = 1.0/domain->yprd;
double invzprd = 1.0/domain->zprd;

32
src/dump_custom.cpp
View File

@ -195,12 +195,12 @@ void DumpCustom::header_binary(int ndump)
{
fwrite(&update->ntimestep,sizeof(int),1,fp);
fwrite(&ndump,sizeof(int),1,fp);
fwrite(&domain->boxxlo,sizeof(double),1,fp);
fwrite(&domain->boxxhi,sizeof(double),1,fp);
fwrite(&domain->boxylo,sizeof(double),1,fp);
fwrite(&domain->boxyhi,sizeof(double),1,fp);
fwrite(&domain->boxzlo,sizeof(double),1,fp);
fwrite(&domain->boxzhi,sizeof(double),1,fp);
fwrite(&boxxlo,sizeof(double),1,fp);
fwrite(&boxxhi,sizeof(double),1,fp);
fwrite(&boxylo,sizeof(double),1,fp);
fwrite(&boxyhi,sizeof(double),1,fp);
fwrite(&boxzlo,sizeof(double),1,fp);
fwrite(&boxzhi,sizeof(double),1,fp);
fwrite(&size_one,sizeof(int),1,fp);
if (multiproc) {
int one = 1;
@ -217,9 +217,9 @@ void DumpCustom::header_item(int ndump)
fprintf(fp,"ITEM: NUMBER OF ATOMS\n");
fprintf(fp,"%d\n",ndump);
fprintf(fp,"ITEM: BOX BOUNDS\n");
fprintf(fp,"%g %g\n",domain->boxxlo,domain->boxxhi);
fprintf(fp,"%g %g\n",domain->boxylo,domain->boxyhi);
fprintf(fp,"%g %g\n",domain->boxzlo,domain->boxzhi);
fprintf(fp,"%g %g\n",boxxlo,boxxhi);
fprintf(fp,"%g %g\n",boxylo,boxyhi);
fprintf(fp,"%g %g\n",boxzlo,boxzhi);
fprintf(fp,"ITEM: ATOMS\n");
}
@ -550,6 +550,7 @@ void DumpCustom::parse_fields(int narg, char **arg)
int i;
for (int iarg = 5; iarg < narg; iarg++) {
i = iarg-5;
if (strcmp(arg[iarg],"tag") == 0) {
pack_choice[i] = &DumpCustom::pack_tag;
vtype[i] = INT;
@ -572,12 +573,18 @@ void DumpCustom::parse_fields(int narg, char **arg)
pack_choice[i] = &DumpCustom::pack_z;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"xs") == 0) {
if (domain->triclinic)
error->all("Cannot dump scaled coords with triclinic box");
pack_choice[i] = &DumpCustom::pack_xs;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"ys") == 0) {
if (domain->triclinic)
error->all("Cannot dump scaled coords with triclinic box");
pack_choice[i] = &DumpCustom::pack_ys;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"zs") == 0) {
if (domain->triclinic)
error->all("Cannot dump scaled coords with triclinic box");
pack_choice[i] = &DumpCustom::pack_zs;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"xu") == 0) {
@ -844,6 +851,13 @@ int DumpCustom::modify_param(int narg, char **arg)
memory->srealloc(thresh_value,(nthresh+1)*sizeof(double),
"dump:thresh_value");
// error check for triclinic box
if (domain->triclinic &&
(strcmp(arg[1],"xs") == 0 || strcmp(arg[1],"ys") == 0 ||
strcmp(arg[1],"zs") == 0))
error->all("Cannot dump scaled coords with triclinic box");
// set keyword type of threshhold
// customize by adding to if statement

23
src/dump_dcd.cpp
View File

@ -15,6 +15,7 @@
Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
------------------------------------------------------------------------- */
#include "math.h"
#include "inttypes.h"
#include "stdio.h"
#include "time.h"
@ -133,14 +134,26 @@ void DumpDCD::write_header(int n)
nframes = 0;
}
// dim[134] = angle cosines of periodic box
// dim[] = size and angle cosines of orthogonal or triclinic box
// dim[1] = alpha = cosine of angle between b and c
// dim[3] = beta = cosine of angle between a and c
// dim[4] = gamma = cosine of angle between a and b
// 48 = 6 doubles
double dim[6];
dim[0] = domain->boxxhi - domain->boxxlo;
dim[2] = domain->boxyhi - domain->boxylo;
dim[5] = domain->boxzhi - domain->boxzlo;
dim[1] = dim[3] = dim[4] = 0.0;
dim[0] = domain->xprd;
dim[2] = domain->yprd;
dim[5] = domain->zprd;
if (domain->triclinic == 0) dim[1] = dim[3] = dim[4] = 0.0;
else {
double *h = domain->h;
double alen = h[0];
double blen = sqrt(h[5]*h[5] + h[1]*h[1]);
double clen = sqrt(h[4]*h[4] + h[3]*h[3] + h[2]*h[2]);
dim[1] = (h[5]*h[4] + h[1]*h[3]) / blen/clen;
dim[3] = (h[0]*h[4]) / alen/clen;
dim[4] = (h[0]*h[5]) / alen/blen;
}
uint32_t out_integer = 48;
fwrite(&out_integer,sizeof(uint32_t),1,fp);

1
src/fix.cpp
View File

@ -36,6 +36,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
restart_global = 0;
restart_peratom = 0;
force_reneighbor = 0;
box_change = 0;
thermo_energy = 0;
pressure_every = 0;

3
src/fix.h
View File

@ -26,6 +26,7 @@ class Fix : protected Pointers {
int restart_global; // 1 if Fix saves global state, 0 if not
int restart_peratom; // 1 if Fix saves peratom state, 0 if not
int force_reneighbor; // 1 if Fix forces reneighboring, 0 if not
int box_change; // 1 if Fix changes box size, 0 if not
int next_reneighbor; // next timestep to force a reneighboring
int thermo_energy; // 1 if ThEng enabled via fix_modify, 0 if not
int nevery; // how often to call an end_of_step fix
@ -81,7 +82,7 @@ class Fix : protected Pointers {
virtual int min_dof() {return 0;}
virtual void min_xinitial(double *) {}
virtual int pack_comm(int, int *, double *, int *) {return 0;}
virtual int pack_comm(int, int *, double *, int, double *) {return 0;}
virtual void unpack_comm(int, int, double *) {}
virtual int pack_reverse_comm(int, int, double *) {return 0;}
virtual void unpack_reverse_comm(int, int *, double *) {}

9
src/fix_ave_spatial.cpp
View File

@ -147,6 +147,15 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
if (delta <= 0.0) error->all("Illegal fix ave/spatial command");
invdelta = 1.0/delta;
if (domain->triclinic) {
if (dim == 0 && (domain->xy != 0.0 || domain->xz != 0.0))
error->all("Cannot (yet) use fix ave/spatial with triclinic box");
if (dim == 1 && (domain->xy != 0.0 || domain->yz != 0.0))
error->all("Cannot (yet) use fix ave/spatial with triclinic box");
if (dim == 2 && (domain->xz != 0.0 || domain->yz != 0.0))
error->all("Cannot (yet) use fix ave/spatial with triclinic box");
}
nvalues = 1;
if (which == COMPUTE) {
int icompute = modify->find_compute(id_compute);

73
src/fix_deposit.cpp
View File

@ -150,7 +150,8 @@ int FixDeposit::setmask()
void FixDeposit::pre_exchange()
{
int flag,flagall;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double coord[3],lamda[3],delx,dely,delz,rsq;
double *newcoord;
// just return if should not be called on this timestep
@ -161,6 +162,15 @@ void FixDeposit::pre_exchange()
double offset = 0.0;
if (rateflag) offset = (update->ntimestep - nfirst) * update->dt * rate;
double *sublo,*subhi;
if (domain->triclinic == 0) {
sublo = domain->sublo;
subhi = domain->subhi;
} else {
sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda;
}
// attempt an insertion until successful
int success = 0;
@ -170,19 +180,19 @@ void FixDeposit::pre_exchange()
// choose random position for new atom within region
xtmp = xlo + random->uniform() * (xhi-xlo);
ytmp = ylo + random->uniform() * (yhi-ylo);
ztmp = zlo + random->uniform() * (zhi-zlo);
while (domain->regions[iregion]->match(xtmp,ytmp,ztmp) == 0) {
xtmp = xlo + random->uniform() * (xhi-xlo);
ytmp = ylo + random->uniform() * (yhi-ylo);
ztmp = zlo + random->uniform() * (zhi-zlo);
coord[0] = xlo + random->uniform() * (xhi-xlo);
coord[1] = ylo + random->uniform() * (yhi-ylo);
coord[2] = zlo + random->uniform() * (zhi-zlo);
while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) {
coord[0] = xlo + random->uniform() * (xhi-xlo);
coord[1] = ylo + random->uniform() * (yhi-ylo);
coord[2] = zlo + random->uniform() * (zhi-zlo);
}
// adjust vertical coord by offset
if (force->dimension == 2) ytmp += offset;
else ztmp += offset;
if (force->dimension == 2) coord[1] += offset;
else coord[2] += offset;
// if global, reset vertical coord to be lo-hi above highest atom
// if local, reset vertical coord to be lo-hi above highest "nearby" atom
@ -204,10 +214,10 @@ void FixDeposit::pre_exchange()
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (localflag) {
delx = xtmp - x[i][0];
dely = ytmp - x[i][1];
delx = coord[0] - x[i][0];
dely = coord[1] - x[i][1];
delz = 0.0;
domain->minimum_image(&delx,&dely,&delz);
domain->minimum_image(delx,dely,delz);
if (force->dimension == 2) rsq = delx*delx;
else rsq = delx*delx + dely*dely;
if (rsq > deltasq) continue;
@ -217,12 +227,12 @@ void FixDeposit::pre_exchange()
MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world);
if (force->dimension == 2)
ytmp = maxall + lo + random->uniform()*(hi-lo);
coord[1] = maxall + lo + random->uniform()*(hi-lo);
else
ztmp = maxall + lo + random->uniform()*(hi-lo);
coord[2] = maxall + lo + random->uniform()*(hi-lo);
}
// now have final xtmp,ytmp,ztmp
// now have final coord
// if distance to any atom is less than near, try again
double **x = atom->x;
@ -230,10 +240,10 @@ void FixDeposit::pre_exchange()
flag = 0;
for (int i = 0; i < nlocal; i++) {
delx = xtmp - x[i][0];
dely = ytmp - x[i][1];
delz = ztmp - x[i][2];
domain->minimum_image(&delx,&dely,&delz);
delx = coord[0] - x[i][0];
dely = coord[1] - x[i][1];
delz = coord[2] - x[i][2];
domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < nearsq) flag = 1;
}
@ -251,21 +261,26 @@ void FixDeposit::pre_exchange()
// if so, add to my list via create_atom()
// initialize info about the atoms
// set group mask to "all" plus fix group
if (domain->triclinic) {
domain->x2lamda(coord,lamda);
newcoord = lamda;
} else newcoord = coord;
flag = 0;
if (xtmp >= domain->subxlo && xtmp < domain->subxhi &&
ytmp >= domain->subylo && ytmp < domain->subyhi &&
ztmp >= domain->subzlo && ztmp < domain->subzhi) flag = 1;
else if (force->dimension == 3 && ztmp >= domain->boxzhi &&
if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] &&
newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1;
else if (force->dimension == 3 && newcoord[2] >= domain->boxzhi &&
comm->myloc[2] == comm->procgrid[2]-1 &&
xtmp >= domain->subxlo && xtmp < domain->subxhi &&
ytmp >= domain->subylo && ytmp < domain->subyhi) flag = 1;
else if (force->dimension == 2 && ytmp >= domain->boxyhi &&
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1;
else if (force->dimension == 2 && newcoord[1] >= domain->boxyhi &&
comm->myloc[1] == comm->procgrid[1]-1 &&
xtmp >= domain->subxlo && xtmp < domain->subxhi) flag = 1;
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1;
if (flag) {
atom->avec->create_atom(ntype,xtmp,ytmp,ztmp,0);
atom->avec->create_atom(ntype,coord,0);
int m = atom->nlocal - 1;
atom->type[m] = ntype;
atom->mask[m] = 1 | groupbit;

2
src/fix_drag.cpp
View File

@ -96,7 +96,7 @@ void FixDrag::post_force(int vflag)
if (!xflag) dx = 0.0;
if (!yflag) dy = 0.0;
if (!zflag) dz = 0.0;
domain->minimum_image(&dx,&dy,&dz);
domain->minimum_image(dx,dy,dz);
r = sqrt(dx*dx + dy*dy + dz*dz);
if (r > delta) {
prefactor = f_mag/r;

10
src/fix_nph.cpp
View File

@ -45,6 +45,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
restart_global = 1;
pressure_every = 1;
box_change = 1;
double p_period[3];
if (strcmp(arg[3],"xyz") == 0) {
@ -119,14 +120,15 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
} else error->all("Illegal fix nph command");
}
// check for periodicity in controlled dimensions
// error checks
if (domain->triclinic) error->all("Cannot use fix nph with triclinic box");
if (p_flag[0] && domain->xperiodic == 0)
error->all("Cannot fix nph on a non-periodic dimension");
error->all("Cannot use fix nph on a non-periodic dimension");
if (p_flag[1] && domain->yperiodic == 0)
error->all("Cannot fix nph on a non-periodic dimension");
error->all("Cannot use fix nph on a non-periodic dimension");
if (p_flag[2] && domain->zperiodic == 0)
error->all("Cannot fix nph on a non-periodic dimension");
error->all("Cannot use fix nph on a non-periodic dimension");
// convert input periods to frequencies

10
src/fix_npt.cpp
View File

@ -45,6 +45,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
restart_global = 1;
pressure_every = 1;
box_change = 1;
t_start = atof(arg[3]);
t_stop = atof(arg[4]);
@ -126,14 +127,15 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
} else error->all("Illegal fix npt command");
}
// check for periodicity in controlled dimensions
// error checks
if (domain->triclinic) error->all("Cannot use fix npt with triclinic box");
if (p_flag[0] && domain->xperiodic == 0)
error->all("Cannot fix npt on a non-periodic dimension");
error->all("Cannot use fix npt on a non-periodic dimension");
if (p_flag[1] && domain->yperiodic == 0)
error->all("Cannot fix npt on a non-periodic dimension");
error->all("Cannot use fix npt on a non-periodic dimension");
if (p_flag[2] && domain->zperiodic == 0)
error->all("Cannot fix npt on a non-periodic dimension");
error->all("Cannot use fix npt on a non-periodic dimension");
// convert input periods to frequencies

3
src/fix_orient_fcc.cpp
View File

@ -419,7 +419,8 @@ double FixOrientFCC::thermo(int n)
/* ---------------------------------------------------------------------- */
int FixOrientFCC::pack_comm(int n, int *list, double *buf, int *pbc_flags)
int FixOrientFCC::pack_comm(int n, int *list, double *buf,
int pbc_flag, double *pair_dist)
{
int i,j,k,id,num;
int *tag = atom->tag;

2
src/fix_orient_fcc.h
View File

@ -45,7 +45,7 @@ class FixOrientFCC : public Fix {
void post_force(int);
void post_force_respa(int, int, int);
double thermo(int);
int pack_comm(int, int *, double *, int *);
int pack_comm(int, int *, double *, int, double *);
void unpack_comm(int, int, double *);
int memory_usage();

21
src/fix_recenter.cpp
View File

@ -135,22 +135,31 @@ void FixRecenter::init()
void FixRecenter::initial_integrate()
{
// target COM
// bounding box around domain works for both orthog and triclinic
double xtarget,ytarget,ztarget;
double *bboxlo,*bboxhi;
if (scaleflag == 2) {
if (domain->triclinic == 0) {
bboxlo = domain->boxlo;
bboxhi = domain->boxhi;
} else {
bboxlo = domain->boxlo_bound;
bboxhi = domain->boxhi_bound;
}
}
if (xinitflag) xtarget = xinit;
else if (scaleflag == 2)
xtarget = domain->boxxlo + xcom*(domain->boxxhi - domain->boxxlo);
else if (scaleflag == 2) xtarget = bboxlo[0] + xcom*(bboxhi[0] - bboxlo[0]);
else xtarget = xcom;
if (yinitflag) ytarget = yinit;
else if (scaleflag == 2)
ytarget = domain->boxylo + ycom*(domain->boxyhi - domain->boxylo);
else if (scaleflag == 2) ytarget = bboxlo[1] + ycom*(bboxhi[1] - bboxlo[1]);
else ytarget = ycom;
if (zinitflag) ztarget = zinit;
else if (scaleflag == 2)
ztarget = domain->boxzlo + zcom*(domain->boxzhi - domain->boxzlo);
else if (scaleflag == 2) ztarget = bboxlo[2] + zcom*(bboxhi[2] - bboxlo[2]);
else ztarget = zcom;
// current COM

19
src/fix_shake.cpp
View File

@ -1911,7 +1911,7 @@ void FixShake::stats()
delx = x[iatom][0] - x[jatom][0];
dely = x[iatom][1] - x[jatom][1];
delz = x[iatom][2] - x[jatom][2];
domain->minimum_image(&delx,&dely,&delz);
domain->minimum_image(delx,dely,delz);
r = sqrt(delx*delx + dely*dely + delz*delz);
m = shake_type[i][j-1];
@ -1931,19 +1931,19 @@ void FixShake::stats()
delx = x[iatom][0] - x[jatom][0];
dely = x[iatom][1] - x[jatom][1];
delz = x[iatom][2] - x[jatom][2];
domain->minimum_image(&delx,&dely,&delz);
domain->minimum_image(delx,dely,delz);
r1 = sqrt(delx*delx + dely*dely + delz*delz);
delx = x[iatom][0] - x[katom][0];
dely = x[iatom][1] - x[katom][1];
delz = x[iatom][2] - x[katom][2];
domain->minimum_image(&delx,&dely,&delz);
domain->minimum_image(delx,dely,delz);
r2 = sqrt(delx*delx + dely*dely + delz*delz);
delx = x[jatom][0] - x[katom][0];
dely = x[jatom][1] - x[katom][1];
delz = x[jatom][2] - x[katom][2];
domain->minimum_image(&delx,&dely,&delz);
domain->minimum_image(delx,dely,delz);
r3 = sqrt(delx*delx + dely*dely + delz*delz);
angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2));
@ -2252,12 +2252,13 @@ void FixShake::post_force_respa(int vflag_in, int ilevel, int iloop)
/* ---------------------------------------------------------------------- */
int FixShake::pack_comm(int n, int *list, double *buf, int *pbc_flags)
int FixShake::pack_comm(int n, int *list, double *buf,
int pbc_flag, double *pbc_dist)
{
int i,j,m;
m = 0;
if (pbc_flags[0] == 0) {
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = xshake[j][0];
@ -2267,9 +2268,9 @@ int FixShake::pack_comm(int n, int *list, double *buf, int *pbc_flags)
} else {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = xshake[j][0] + pbc_flags[1]*domain->xprd;
buf[m++] = xshake[j][1] + pbc_flags[2]*domain->yprd;
buf[m++] = xshake[j][2] + pbc_flags[3]*domain->zprd;
buf[m++] = xshake[j][0] + pbc_dist[0];
buf[m++] = xshake[j][1] + pbc_dist[1];
buf[m++] = xshake[j][2] + pbc_dist[2];
}
}
return 3;

2
src/fix_shake.h
View File

@ -34,7 +34,7 @@ class FixShake : public Fix {
void copy_arrays(int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
int pack_comm(int, int *, double *, int *);
int pack_comm(int, int *, double *, int, double *);
void unpack_comm(int, int, double *);
int dof(int);

6
src/fix_uniaxial.cpp
View File

@ -47,8 +47,11 @@ FixUniaxial::FixUniaxial(LAMMPS *lmp, int narg, char **arg) :
lambda_final = atof(arg[5]);
if (lambda_final <= 0) error->all("Illegal fix uniaxial command");
if (domain->nonperiodic)
error->all("Cannot fix uniaxial on non-periodic system");
error->all("Cannot use fix uniaxial on non-periodic system");
if (domain->triclinic)
error->all("Cannot use fix uniaxial with triclinic box");
nrigid = 0;
rfix = NULL;
@ -204,4 +207,3 @@ void FixUniaxial::end_of_step()
if (kspace_flag) force->kspace->setup();
}

6
src/fix_volume_rescale.cpp
View File

@ -31,6 +31,9 @@ FixVolRescale::FixVolRescale(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all("Illegal fix volume/rescale command");
box_change = 1;
nevery = atoi(arg[3]);
if (nevery <= 0) error->all("Illegal fix volume/rescale command");
@ -64,6 +67,9 @@ FixVolRescale::FixVolRescale(LAMMPS *lmp, int narg, char **arg) :
} else error->all("Illegal fix volume/rescale command");
}
if (domain->triclinic)
error->all("Cannot use fix volume/rescale with triclinic box");
nrigid = 0;
rfix = NULL;
}

8
src/fix_wall_lj126.cpp
View File

@ -20,6 +20,7 @@
#include "string.h"
#include "fix_wall_lj126.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
#include "output.h"
#include "respa.h"
@ -67,6 +68,13 @@ FixWallLJ126::FixWallLJ126(LAMMPS *lmp, int narg, char **arg) :
double r2inv = 1.0/(cutoff*cutoff);
double r6inv = r2inv*r2inv*r2inv;
offset = r6inv*(coeff3*r6inv - coeff4);
if (dim == 0 && domain->xperiodic)
error->all("Cannot use wall in periodic dimension");
if (dim == 1 && domain->yperiodic)
error->all("Cannot use wall in periodic dimension");
if (dim == 2 && domain->zperiodic)
error->all("Cannot use wall in periodic dimension");
}
/* ---------------------------------------------------------------------- */

8
src/fix_wall_lj93.cpp
View File

@ -16,6 +16,7 @@
#include "string.h"
#include "fix_wall_lj93.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
#include "output.h"
#include "respa.h"
@ -64,6 +65,13 @@ FixWallLJ93::FixWallLJ93(LAMMPS *lmp, int narg, char **arg) :
double r2inv = rinv*rinv;
double r4inv = r2inv*r2inv;
offset = coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv;
if (dim == 0 && domain->xperiodic)
error->all("Cannot use wall in periodic dimension");
if (dim == 1 && domain->yperiodic)
error->all("Cannot use wall in periodic dimension");
if (dim == 2 && domain->zperiodic)
error->all("Cannot use wall in periodic dimension");
}
/* ---------------------------------------------------------------------- */

7
src/fix_wall_reflect.cpp
View File

@ -38,6 +38,13 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg],"zhi") == 0) zhiflag = 1;
else error->all("Illegal fix wall/reflect command");
}
if ((xloflag || xhiflag) && domain->xperiodic)
error->all("Cannot use wall in periodic dimension");
if ((yloflag || yhiflag) && domain->yperiodic)
error->all("Cannot use wall in periodic dimension");
if ((zloflag || zhiflag) && domain->zperiodic)
error->all("Cannot use wall in periodic dimension");
}
/* ---------------------------------------------------------------------- */

76
src/lattice.cpp
View File

@ -27,7 +27,7 @@ using namespace LAMMPS_NS;
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
#define BIG 1.0e30
enum{NONE,SC,BCC,FCC,DIAMOND,SQ,SQ2,HEX,USER};
enum{NONE,SC,BCC,FCC,DIAMOND,SQ,SQ2,HEX,CUSTOM};
/* ---------------------------------------------------------------------- */
@ -45,7 +45,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
else if (strcmp(arg[0],"sq") == 0) style = SQ;
else if (strcmp(arg[0],"sq2") == 0) style = SQ2;
else if (strcmp(arg[0],"hex") == 0) style = HEX;
else if (strcmp(arg[0],"user") == 0) style = USER;
else if (strcmp(arg[0],"custom") == 0) style = CUSTOM;
else error->all("Illegal lattice command");
if (style == NONE) {
@ -54,7 +54,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
}
// check that lattice matches dimension
// style USER can be either 2d or 3d
// style CUSTOM can be either 2d or 3d
int dimension = force->dimension;
if (dimension == 2) {
@ -74,7 +74,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
// set basis atoms for each style
// x,y,z = fractional coords within unit cell
// style USER will be defined by optional args
// style CUSTOM will be defined by optional args
nbasis = 0;
basis = NULL;
@ -116,6 +116,8 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
orienty[0] = 0; orienty[1] = 1; orienty[2] = 0;
orientz[0] = 0; orientz[1] = 0; orientz[2] = 1;
int spaceflag = 0;
a1[0] = 1.0; a1[1] = 0.0; a1[2] = 0.0;
a2[0] = 0.0; a2[1] = 1.0; a2[2] = 0.0;
a3[0] = 0.0; a3[1] = 0.0; a3[2] = 1.0;
@ -153,26 +155,34 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
orient[2] = atoi(arg[iarg+4]);
iarg += 5;
} else if (strcmp(arg[iarg],"spacings") == 0) {
if (iarg+4 > narg) error->all("Illegal lattice command");
spaceflag = 1;
xlattice = atoi(arg[iarg+1]);
ylattice = atoi(arg[iarg+2]);
zlattice = atoi(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"a1") == 0) {
if (iarg+4 > narg) error->all("Illegal lattice command");
if (style != USER)
error->all("Invalid option in lattice command for non-user style");
if (style != CUSTOM)
error->all("Invalid option in lattice command for non-custom style");
a1[0] = atof(arg[iarg+1]);
a1[1] = atof(arg[iarg+2]);
a1[2] = atof(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"a2") == 0) {
if (iarg+4 > narg) error->all("Illegal lattice command");
if (style != USER)
error->all("Invalid option in lattice command for non-user style");
if (style != CUSTOM)
error->all("Invalid option in lattice command for non-custom style");
a2[0] = atof(arg[iarg+1]);
a2[1] = atof(arg[iarg+2]);
a2[2] = atof(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"a3") == 0) {
if (iarg+4 > narg) error->all("Illegal lattice command");
if (style != USER)
error->all("Invalid option in lattice command for non-user style");
if (style != CUSTOM)
error->all("Invalid option in lattice command for non-custom style");
a3[0] = atof(arg[iarg+1]);
a3[1] = atof(arg[iarg+2]);
a3[2] = atof(arg[iarg+3]);
@ -180,8 +190,8 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
} else if (strcmp(arg[iarg],"basis") == 0) {
if (iarg+4 > narg) error->all("Illegal lattice command");
if (style != USER)
error->all("Invalid option in lattice command for non-user style");
if (style != CUSTOM)
error->all("Invalid option in lattice command for non-custom style");
double x = atof(arg[iarg+1]);
double y = atof(arg[iarg+2]);
double z = atof(arg[iarg+3]);
@ -212,6 +222,11 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
error->all("Lattice settings are not compatible with 2d simulation");
}
if (spaceflag) {
if (xlattice <= 0.0 || ylattice <= 0.0 || zlattice <= 0.0)
error->all("Lattice spacings are invalid");
}
// reset scale for LJ units (input scale is rho*)
// scale = (Nbasis/(Vprimitive * rho*)) ^ (1/dim)
@ -232,23 +247,30 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
// set xlattice,ylattice,zlattice to 0.0 initially
// since bbox uses them to shift origin (irrelevant for this computation)
double xmin,ymin,zmin,xmax,ymax,zmax;
xmin = ymin = zmin = BIG;
xmax = ymax = zmax = -BIG;
xlattice = ylattice = zlattice = 0.0;
if (spaceflag == 0) {
double xmin,ymin,zmin,xmax,ymax,zmax;
xmin = ymin = zmin = BIG;
xmax = ymax = zmax = -BIG;
xlattice = ylattice = zlattice = 0.0;
bbox(0,0.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,1.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,0.0,1.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,1.0,1.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,0.0,0.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,1.0,0.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,0.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,1.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,0.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,1.0,0.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,0.0,1.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,1.0,1.0,0.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,0.0,0.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,1.0,0.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,0.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax);
bbox(0,1.0,1.0,1.0,xmin,ymin,zmin,xmax,ymax,zmax);
xlattice = xmax - xmin;
ylattice = ymax - ymin;
zlattice = zmax - zmin;
xlattice = xmax - xmin;
ylattice = ymax - ymin;
zlattice = zmax - zmin;
} else {
xlattice *= scale;
ylattice *= scale;
zlattice *= scale;
}
// print lattice spacings

15
src/neigh_full.cpp
View File

@ -131,16 +131,9 @@ void Neighbor::full_bin()
// only skip i = j
for (k = 0; k < nstencil_full; k++) {
j = binhead[ibin+stencil_full[k]];
while (j >= 0) {
if (i == j) {
j = bins[j];
continue;
}
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
}
for (j = binhead[ibin+stencil_full[k]]; j >= 0; j = bins[j]) {
if (i == j) continue;
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
jtype = type[j];
delx = xtmp - x[j][0];
@ -154,8 +147,6 @@ void Neighbor::full_bin()
if (which == 0) neighptr[n++] = j;
else if (which > 0) neighptr[n++] = which*nall + j;
}
j = bins[j];
}
}

127
src/neigh_gran.cpp
View File

@ -290,17 +290,9 @@ void Neighbor::granular_bin_no_newton()
// stores own/ghost pairs on both procs
for (k = 0; k < nstencil; k++) {
j = binhead[ibin+stencil[k]];
while (j >= 0) {
if (j <= i) {
j = bins[j];
continue;
}
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
}
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (j <= i) continue;
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
@ -337,8 +329,6 @@ void Neighbor::granular_bin_no_newton()
n++;
}
j = bins[j];
}
}
@ -358,8 +348,8 @@ void Neighbor::granular_bin_no_newton()
granular particles
binned neighbor list construction with full Newton's 3rd law
no shear history is allowed for this option
every pair stored exactly once by some processor
each owned atom i checks its own bin and other bins in Newton stencil
every pair stored exactly once by some processor
------------------------------------------------------------------------- */
void Neighbor::granular_bin_newton()
@ -405,19 +395,13 @@ void Neighbor::granular_bin_newton()
// if j is owned atom, store it, since j is beyond i in linked list
// if j is ghost, only store if j coords are "above and to the right" of i
j = bins[i];
while (j >= 0) {
if (j >= nlocal) {
if ((x[j][2] < ztmp) || (x[j][2] == ztmp && x[j][1] < ytmp) ||
(x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)) {
j = bins[j];
continue;
}
}
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
for (j = bins[i]; j >= 0; j = bins[j]) {
if (j >= nlocal) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
}
delx = xtmp - x[j][0];
@ -428,20 +412,93 @@ void Neighbor::granular_bin_newton()
cutsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutsq) neighptr[n++] = j;
j = bins[j];
}
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
for (k = 0; k < nstencil; k++) {
j = binhead[ibin+stencil[k]];
while (j >= 0) {
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
}
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutsq) neighptr[n++] = j;
}
}
firstneigh[i] = neighptr;
numneigh[i] = n;
npnt += n;
if (npnt >= pgsize)
error->one("Neighbor list overflow, boost neigh_modify one or page");
}
}
/* ----------------------------------------------------------------------
granular particles
binned neighbor list construction with Newton's 3rd law for triclinic
no shear history is allowed for this option
each owned atom i checks its own bin and other bins in triclinic stencil
every pair stored exactly once by some processor
------------------------------------------------------------------------- */
void Neighbor::granular_bin_newton_tri()
{
int i,j,k,n,ibin;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutsq;
int *neighptr;
// bin local & ghost atoms
bin_atoms();
// loop over each atom, storing neighbors
double **x = atom->x;
double *radius = atom->radius;
int *type = atom->type;
int *mask = atom->mask;
int *molecule = atom->molecule;
int nlocal = atom->nlocal;
int npage = 0;
int npnt = 0;
for (i = 0; i < nlocal; i++) {
if (pgsize - npnt < oneatom) {
npnt = 0;
npage++;
if (npage == maxpage) add_pages(npage);
}
n = 0;
neighptr = &pages[npage][npnt];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
radi = radius[i];
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and <= x)
// this excludes self-self interaction
ibin = coord2bin(x[i]);
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue;
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
@ -451,8 +508,6 @@ void Neighbor::granular_bin_newton()
cutsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutsq) neighptr[n++] = j;
j = bins[j];
}
}

131
src/neigh_half.cpp
View File

@ -219,17 +219,9 @@ void Neighbor::half_bin_no_newton()
// stores own/ghost pairs on both procs
for (k = 0; k < nstencil; k++) {
j = binhead[ibin+stencil[k]];
while (j >= 0) {
if (j <= i) {
j = bins[j];
continue;
}
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
}
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (j <= i) continue;
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
jtype = type[j];
delx = xtmp - x[j][0];
@ -243,8 +235,6 @@ void Neighbor::half_bin_no_newton()
if (which == 0) neighptr[n++] = j;
else if (which > 0) neighptr[n++] = which*nall + j;
}
j = bins[j];
}
}
@ -305,20 +295,14 @@ void Neighbor::half_bin_newton()
// if j is owned atom, store it, since j is beyond i in linked list
// if j is ghost, only store if j coords are "above and to the right" of i
j = bins[i];
while (j >= 0) {
for (j = bins[i]; j >= 0; j = bins[j]) {
if (j >= nlocal) {
if ((x[j][2] < ztmp) || (x[j][2] == ztmp && x[j][1] < ytmp) ||
(x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)) {
j = bins[j];
continue;
}
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
}
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
jtype = type[j];
delx = xtmp - x[j][0];
@ -332,20 +316,14 @@ void Neighbor::half_bin_newton()
if (which == 0) neighptr[n++] = j;
else if (which > 0) neighptr[n++] = which*nall + j;
}
j = bins[j];
}
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
for (k = 0; k < nstencil; k++) {
j = binhead[ibin+stencil[k]];
while (j >= 0) {
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
}
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
jtype = type[j];
delx = xtmp - x[j][0];
@ -359,8 +337,87 @@ void Neighbor::half_bin_newton()
if (which == 0) neighptr[n++] = j;
else if (which > 0) neighptr[n++] = which*nall + j;
}
}
}
j = bins[j];
firstneigh[i] = neighptr;
numneigh[i] = n;
npnt += n;
if (npnt >= pgsize)
error->one("Neighbor list overflow, boost neigh_modify one or page");
}
}
/* ----------------------------------------------------------------------
binned neighbor list construction with Newton's 3rd law for triclinic
each owned atom i checks its own bin and other bins in triclinic stencil
every pair stored exactly once by some processor
------------------------------------------------------------------------- */
void Neighbor::half_bin_newton_tri()
{
int i,j,k,n,itype,jtype,ibin,which;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr;
// bin local & ghost atoms
bin_atoms();
// loop over each atom, storing neighbors
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
int *molecule = atom->molecule;
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
int molecular = atom->molecular;
int npage = 0;
int npnt = 0;
for (i = 0; i < nlocal; i++) {
if (pgsize - npnt < oneatom) {
npnt = 0;
npage++;
if (npage == maxpage) add_pages(npage);
}
neighptr = &pages[npage][npnt];
n = 0;
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and <= x)
// this excludes self-self interaction
ibin = coord2bin(x[i]);
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue;
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) which = find_special(i,j);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (which > 0) neighptr[n++] = which*nall + j;
}
}
}
@ -460,11 +517,9 @@ void Neighbor::half_full_newton()
if (j < nlocal) {
if (i > j) continue;
} else {
if ((x[j][2] < ztmp) || (x[j][2] == ztmp && x[j][1] < ytmp) ||
(x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)) {
j = bins[j];
continue;
}
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
neighptr[n++] = j;
}

177
src/neigh_respa.cpp
View File

@ -342,17 +342,9 @@ void Neighbor::respa_bin_no_newton()
// stores own/ghost pairs on both procs
for (k = 0; k < nstencil; k++) {
j = binhead[ibin+stencil[k]];
while (j >= 0) {
if (j <= i) {
j = bins[j];
continue;
}
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
}
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (j <= i) continue;
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
jtype = type[j];
delx = xtmp - x[j][0];
@ -376,8 +368,6 @@ void Neighbor::respa_bin_no_newton()
else if (which > 0) neighptr_middle[n_middle++] = which*nall + j;
}
}
j = bins[j];
}
}
@ -406,8 +396,8 @@ void Neighbor::respa_bin_no_newton()
/* ----------------------------------------------------------------------
multiple respa lists
binned neighbor list construction with full Newton's 3rd law
every pair stored exactly once by some processor
each owned atom i checks its own bin and other bins in Newton stencil
every pair stored exactly once by some processor
------------------------------------------------------------------------- */
void Neighbor::respa_bin_newton()
@ -477,19 +467,12 @@ void Neighbor::respa_bin_newton()
// if j is owned atom, store it, since j is beyond i in linked list
// if j is ghost, only store if j coords are "above and to the right" of i
j = bins[i];
while (j >= 0) {
for (j = bins[i]; j >= 0; j = bins[j]) {
if (j >= nlocal) {
if ((x[j][2] < ztmp) || (x[j][2] == ztmp && x[j][1] < ytmp) ||
(x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)) {
j = bins[j];
continue;
}
}
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
}
jtype = type[j];
@ -514,20 +497,144 @@ void Neighbor::respa_bin_newton()
else if (which > 0) neighptr_middle[n_middle++] = which*nall + j;
}
}
j = bins[j];
}
// loop over all atoms in other bins in stencil, store every pair
ibin = coord2bin(x[i]);
for (k = 0; k < nstencil; k++) {
j = binhead[ibin+stencil[k]];
while (j >= 0) {
if (exclude && exclusion(i,j,type,mask,molecule)) {
j = bins[j];
continue;
}
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) which = find_special(i,j);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (which > 0) neighptr[n++] = which*nall + j;
if (rsq < cut_inner_sq) {
if (which == 0) neighptr_inner[n_inner++] = j;
else if (which > 0) neighptr_inner[n_inner++] = which*nall + j;
}
if (respa == 2 && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) {
if (which == 0) neighptr_middle[n_middle++] = j;
else if (which > 0) neighptr_middle[n_middle++] = which*nall + j;
}
}
}
}
firstneigh[i] = neighptr;
numneigh[i] = n;
npnt += n;
if (npnt >= pgsize)
error->one("Neighbor list overflow, boost neigh_modify one or page");
firstneigh_inner[i] = neighptr_inner;
numneigh_inner[i] = n_inner;
npnt_inner += n_inner;
if (npnt_inner >= pgsize)
error->one("Neighbor list overflow, boost neigh_modify one or page");
if (respa == 2) {
firstneigh_middle[i] = neighptr_middle;
numneigh_middle[i] = n_middle;
npnt_middle += n_middle;
if (npnt_middle >= pgsize)
error->one("Neighbor list overflow, boost neigh_modify one or page");
}
}
}
/* ----------------------------------------------------------------------
multiple respa lists
binned neighbor list construction with Newton's 3rd law for triclinic
each owned atom i checks its own bin and other bins in triclinic stencil
every pair stored exactly once by some processor
------------------------------------------------------------------------- */
void Neighbor::respa_bin_newton_tri()
{
int i,j,k,itype,jtype,ibin,which;
int n_inner,n_middle,n;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr_inner;
int *neighptr_middle;
int *neighptr;
// bin local & ghost atoms
bin_atoms();
// loop over each atom, storing neighbors
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
int *molecule = atom->molecule;
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
int molecular = atom->molecular;
int npage = 0;
int npnt = 0;
int npage_inner = 0;
int npnt_inner = 0;
int npage_middle = 0;
int npnt_middle = 0;
for (i = 0; i < nlocal; i++) {
if (pgsize - npnt < oneatom) {
npnt = 0;
npage++;
if (npage == maxpage) add_pages(npage);
}
neighptr = &pages[npage][npnt];
n = 0;
if (pgsize - npnt_inner < oneatom) {
npnt_inner = 0;
npage_inner++;
if (npage_inner == maxpage_inner) add_pages_inner(npage_inner);
}
neighptr_inner = &pages_inner[npage_inner][npnt_inner];
n_inner = 0;
if (respa == 2) {
if (pgsize - npnt_middle < oneatom) {
npnt_middle = 0;
npage_middle++;
if (npage_middle == maxpage_middle) add_pages_middle(npage_middle);
}
neighptr_middle = &pages_middle[npage_middle][npnt_middle];
n_middle = 0;
}
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and <= x)
// this excludes self-self interaction
ibin = coord2bin(x[i]);
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] <= xtmp) continue;
if (exclude && exclusion(i,j,type,mask,molecule)) continue;
jtype = type[j];
delx = xtmp - x[j][0];
@ -551,8 +658,6 @@ void Neighbor::respa_bin_newton()
else if (which > 0) neighptr_middle[n_middle++] = which*nall + j;
}
}
j = bins[j];
}
}

127
src/neighbor.cpp
View File

@ -206,6 +206,7 @@ void Neighbor::init()
int i,j,m,n;
ncalls = ndanger = 0;
triclinic = domain->triclinic;
// error check
@ -544,7 +545,7 @@ void Neighbor::init()
cut_middle_inside_sq = (cut_respa[0] - skin) * (cut_respa[0] - skin);
}
// set ptrs to correct half and full build functions
// set ptrs to correct half, full, triclinic build functions
// cannot combine granular and rRESPA
if (half) {
@ -556,6 +557,8 @@ void Neighbor::init()
} else if (style == 1) {
if (force->newton_pair == 0)
half_build = &Neighbor::granular_bin_no_newton;
else if (triclinic)
half_build = &Neighbor::granular_bin_newton_tri;
else half_build = &Neighbor::granular_bin_newton;
}
} else if (respa) {
@ -566,6 +569,8 @@ void Neighbor::init()
} else if (style == 1) {
if (force->newton_pair == 0)
half_build = &Neighbor::respa_bin_no_newton;
else if (triclinic)
half_build = &Neighbor::respa_bin_newton_tri;
else half_build = &Neighbor::respa_bin_newton;
}
} else {
@ -583,6 +588,8 @@ void Neighbor::init()
else half_build = &Neighbor::half_bin_no_newton;
} else {
if (full_every) half_build = &Neighbor::half_full_newton;
else if (triclinic)
half_build = &Neighbor::half_bin_newton_tri;
else half_build = &Neighbor::half_bin_newton;
}
}
@ -869,6 +876,7 @@ void Neighbor::build_full()
binsize = 1/2 of cutoff is roughly optimal
prd must be filled exactly by integer # of bins
so procs on both sides of PBC see same bin boundary
triclinic is an exception, since stencil & neigh list built differently
mbinlo = lowest global bin any of my ghost atoms could fall into
mbinhi = highest global bin any of my ghost atoms could fall into
mbin = number of bins I need in a dimension
@ -878,31 +886,53 @@ void Neighbor::build_full()
void Neighbor::setup_bins()
{
double cutneighinv = 1.0/cutneigh;
// bbox = bounding box of entire domain
// bboxlo,bboxhi = bounding box of proc's sub-domain
// for triclinic, bounding box surrounds all 8 corner pts
double bbox[3];
double *bboxlo,*bboxhi;
if (triclinic == 0) {
boxlo = domain->boxlo;
boxhi = domain->boxhi;
bboxlo = domain->sublo;
bboxhi = domain->subhi;
} else {
boxlo = domain->boxlo_bound;
boxhi = domain->boxhi_bound;
bboxlo = domain->sublo_bound;
bboxhi = domain->subhi_bound;
}
bbox[0] = boxhi[0] - boxlo[0];
bbox[1] = boxhi[1] - boxlo[1];
bbox[2] = boxhi[2] - boxlo[2];
// test for too many global bins in any dimension due to huge domain
if (2.0*domain->xprd*cutneighinv > INT_MAX ||
2.0*domain->yprd*cutneighinv > INT_MAX ||
2.0*domain->zprd*cutneighinv > INT_MAX)
double cutneighinv = 1.0/cutneigh;
if (2.0*bbox[0]*cutneighinv > INT_MAX || 2.0*bbox[1]*cutneighinv > INT_MAX ||
2.0*bbox[2]*cutneighinv > INT_MAX)
error->all("Domain too large for neighbor bins");
// divide box into bins
// optimal size is roughly 1/2 the cutoff
nbinx = static_cast<int> (2.0*domain->xprd*cutneighinv);
nbiny = static_cast<int> (2.0*domain->yprd*cutneighinv);
nbinx = static_cast<int> (2.0*bbox[0]*cutneighinv);
nbiny = static_cast<int> (2.0*bbox[1]*cutneighinv);
if (force->dimension == 3)
nbinz = static_cast<int> (2.0*domain->zprd*cutneighinv);
nbinz = static_cast<int> (2.0*bbox[2]*cutneighinv);
else nbinz = 1;
if (nbinx == 0) nbinx = 1;
if (nbiny == 0) nbiny = 1;
if (nbinz == 0) nbinz = 1;
binsizex = domain->xprd/nbinx;
binsizey = domain->yprd/nbiny;
binsizez = domain->zprd/nbinz;
binsizex = bbox[0]/nbinx;
binsizey = bbox[1]/nbiny;
binsizez = bbox[2]/nbinz;
bininvx = 1.0 / binsizex;
bininvy = 1.0 / binsizey;
@ -915,23 +945,23 @@ void Neighbor::setup_bins()
double coord;
int mbinxhi,mbinyhi,mbinzhi;
coord = domain->subxlo - cutneigh - SMALL*domain->xprd;
mbinxlo = static_cast<int> ((coord-domain->boxxlo)*bininvx);
coord = bboxlo[0] - cutneigh - SMALL*bbox[0];
mbinxlo = static_cast<int> ((coord-boxlo[0])*bininvx);
if (coord < 0.0) mbinxlo = mbinxlo - 1;
coord = domain->subxhi + cutneigh + SMALL*domain->xprd;
mbinxhi = static_cast<int> ((coord-domain->boxxlo)*bininvx);
coord = bboxhi[0] + cutneigh + SMALL*bbox[0];
mbinxhi = static_cast<int> ((coord-boxlo[0])*bininvx);
coord = domain->subylo - cutneigh - SMALL*domain->yprd;
mbinylo = static_cast<int> ((coord-domain->boxylo)*bininvy);
coord = bboxlo[1] - cutneigh - SMALL*bbox[1];
mbinylo = static_cast<int> ((coord-boxlo[1])*bininvy);
if (coord < 0.0) mbinylo = mbinylo - 1;
coord = domain->subyhi + cutneigh + SMALL*domain->yprd;
mbinyhi = static_cast<int> ((coord-domain->boxylo)*bininvy);
coord = bboxhi[1] + cutneigh + SMALL*bbox[1];
mbinyhi = static_cast<int> ((coord-boxlo[1])*bininvy);
coord = domain->subzlo - cutneigh - SMALL*domain->zprd;
mbinzlo = static_cast<int> ((coord-domain->boxzlo)*bininvz);
coord = bboxlo[2] - cutneigh - SMALL*bbox[2];
mbinzlo = static_cast<int> ((coord-boxlo[2])*bininvz);
if (coord < 0.0) mbinzlo = mbinzlo - 1;
coord = domain->subzhi + cutneigh + SMALL*domain->zprd;
mbinzhi = static_cast<int> ((coord-domain->boxzlo)*bininvz);
coord = bboxhi[2] + cutneigh + SMALL*bbox[2];
mbinzhi = static_cast<int> ((coord-boxlo[2])*bininvz);
// extend bins by 1 to insure stencil extent is included
@ -965,10 +995,15 @@ void Neighbor::setup_bins()
// is within neighbor cutoff
// next(xyz) = how far the stencil could possibly extend
// for partial Newton (newton = 0)
// stencil is all surrounding bins including self
// stencil is all surrounding bins
// stencil includes self
// for triclinic
// stencil is all bins in z-plane of self and above, but not below
// in 2d is all bins in y-plane of self and above, but not below
// stencil includes self
// for full Newton (newton = 1)
// stencil is bins to the "upper right" of central bin
// stencil does NOT include self
// stencil does not include self
// for full neighbor list (full = 1)
// stencil is all surrounding bins including self, regardless of Newton
// stored in stencil_full
@ -999,6 +1034,12 @@ void Neighbor::setup_bins()
for (i = -nextx; i <= nextx; i++)
if (bin_distance(i,j,k) < cutsq)
stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
} else if (triclinic) {
for (k = 0; k <= nextz; k++)
for (j = -nexty; j <= nexty; j++)
for (i = -nextx; i <= nextx; i++)
if (bin_distance(i,j,k) < cutsq)
stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
} else {
for (k = 0; k <= nextz; k++)
for (j = -nexty; j <= nexty; j++)
@ -1013,6 +1054,11 @@ void Neighbor::setup_bins()
for (i = -nextx; i <= nextx; i++)
if (bin_distance(i,j,0) < cutsq)
stencil[nstencil++] = j*mbinx + i;
} else if (triclinic) {
for (j = 0; j <= nexty; j++)
for (i = -nextx; i <= nextx; i++)
if (bin_distance(i,j,0) < cutsq)
stencil[nstencil++] = j*mbinx + i;
} else {
for (j = 0; j <= nexty; j++)
for (i = -nextx; i <= nextx; i++)
@ -1373,32 +1419,33 @@ void Neighbor::bin_atoms()
take special care to insure ghosts are put in correct bins
this is necessary so that both procs on either side of PBC
treat a pair of atoms straddling the PBC in a consistent way
triclinic is an exception, since stencil & neigh list built differently
------------------------------------------------------------------------- */
int Neighbor::coord2bin(double *x)
{
int ix,iy,iz;
if (x[0] >= domain->boxxhi)
ix = static_cast<int> ((x[0]-domain->boxxhi)*bininvx) + nbinx - mbinxlo;
else if (x[0] >= domain->boxxlo)
ix = static_cast<int> ((x[0]-domain->boxxlo)*bininvx) - mbinxlo;
if (x[0] >= boxhi[0])
ix = static_cast<int> ((x[0]-boxhi[0])*bininvx) + nbinx - mbinxlo;
else if (x[0] >= boxlo[0])
ix = static_cast<int> ((x[0]-boxlo[0])*bininvx) - mbinxlo;
else
ix = static_cast<int> ((x[0]-domain->boxxlo)*bininvx) - mbinxlo - 1;
ix = static_cast<int> ((x[0]-boxlo[0])*bininvx) - mbinxlo - 1;
if (x[1] >= domain->boxyhi)
iy = static_cast<int> ((x[1]-domain->boxyhi)*bininvy) + nbiny - mbinylo;
else if (x[1] >= domain->boxylo)
iy = static_cast<int> ((x[1]-domain->boxylo)*bininvy) - mbinylo;
if (x[1] >= boxhi[1])
iy = static_cast<int> ((x[1]-boxhi[1])*bininvy) + nbiny - mbinylo;
else if (x[1] >= boxlo[1])
iy = static_cast<int> ((x[1]-boxlo[1])*bininvy) - mbinylo;
else
iy = static_cast<int> ((x[1]-domain->boxylo)*bininvy) - mbinylo - 1;
iy = static_cast<int> ((x[1]-boxlo[1])*bininvy) - mbinylo - 1;
if (x[2] >= domain->boxzhi)
iz = static_cast<int> ((x[2]-domain->boxzhi)*bininvz) + nbinz - mbinzlo;
else if (x[2] >= domain->boxzlo)
iz = static_cast<int> ((x[2]-domain->boxzlo)*bininvz) - mbinzlo;
if (x[2] >= boxhi[2])
iz = static_cast<int> ((x[2]-boxhi[2])*bininvz) + nbinz - mbinzlo;
else if (x[2] >= boxlo[2])
iz = static_cast<int> ((x[2]-boxlo[2])*bininvz) - mbinzlo;
else
iz = static_cast<int> ((x[2]-domain->boxzlo)*bininvz) - mbinzlo - 1;
iz = static_cast<int> ((x[2]-boxlo[2])*bininvz) - mbinzlo - 1;
return (iz*mbiny*mbinx + iy*mbinx + ix + 1);
}

9
src/neighbor.h
View File

@ -117,6 +117,9 @@ class Neighbor : protected Pointers {
double binsizex,binsizey,binsizez; // bin sizes and inverse sizes
double bininvx,bininvy,bininvz;
int triclinic; // 0 if domain is orthog, 1 if triclinic
double *boxlo,*boxhi; // copies of domain bounds
int **pages; // half neighbor list pages
int maxpage; // # of half pages currently allocated
@ -217,6 +220,12 @@ class Neighbor : protected Pointers {
int coord2bin(double *); // mapping atom coord to a bin
int find_special(int, int); // look for special neighbor
int exclusion(int, int, int *, int *, int *); // test for pair exclusion
// PJV
void granular_bin_newton_tri();
void respa_bin_newton_tri();
void half_bin_newton_tri();
};
}

2
src/pair.h
View File

@ -79,7 +79,7 @@ class Pair : protected Pointers {
virtual void write_restart_settings(FILE *) {}
virtual void read_restart_settings(FILE *) {}
virtual int pack_comm(int, int *, double *, int *) {return 0;}
virtual int pack_comm(int, int *, double *, int, double *) {return 0;}
virtual void unpack_comm(int, int, double *) {}
virtual int pack_reverse_comm(int, int, double *) {return 0;}
virtual void unpack_reverse_comm(int, int *, double *) {}

11
src/read_data.cpp
View File

@ -11,6 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "mpi.h"
#include "string.h"
#include "stdlib.h"
@ -114,6 +115,7 @@ void ReadData::command(int narg, char **arg)
atom->allocate_type_arrays();
atom->avec->grow(n);
domain->print_box(" ");
domain->set_initial_box();
domain->set_global_box();
comm->set_procs();
@ -356,7 +358,10 @@ void ReadData::header(int flag)
sscanf(line,"%lg %lg",&domain->boxylo,&domain->boxyhi);
else if (strstr(line,"zlo zhi"))
sscanf(line,"%lg %lg",&domain->boxzlo,&domain->boxzhi);
else break;
else if (strstr(line,"xy xz yz")) {
domain->triclinic = 1;
sscanf(line,"%lg %lg %lg",&domain->xy,&domain->xz,&domain->yz);
} else break;
}
// error check on consistency of header values
@ -382,10 +387,6 @@ void ReadData::header(int flag)
error->one("Dihedrals defined but no dihedral types");
if (atom->nimpropers > 0 && atom->nimpropertypes <= 0)
error->one("Impropers defined but no improper types");
if (domain->boxxlo >= domain->boxxhi || domain->boxylo >= domain->boxyhi ||
domain->boxzlo >= domain->boxzhi)
error->one("Box bounds are invalid");
}
/* ----------------------------------------------------------------------

44
src/read_restart.cpp
View File

@ -101,6 +101,7 @@ void ReadRestart::command(int narg, char **arg)
atom->allocate_type_arrays();
atom->avec->grow(n);
domain->print_box(" ");
domain->set_initial_box();
domain->set_global_box();
comm->set_procs();
@ -124,21 +125,25 @@ void ReadRestart::command(int narg, char **arg)
// proc 0 reads chunks from series of files (nprocs_file)
// proc 0 bcasts each chunk to other procs
// each proc unpacks the atoms, saving ones in it's sub-domain
// check for atom in sub-domain is different for orthogonal vs triclinic box
AtomVec *avec = atom->avec;
int triclinic = domain->triclinic;
int maxbuf = 0;
double *buf = NULL;
double subxlo = domain->subxlo;
double subxhi = domain->subxhi;
double subylo = domain->subylo;
double subyhi = domain->subyhi;
double subzlo = domain->subzlo;
double subzhi = domain->subzhi;
double *x,lamda[3];
double *coord,*sublo,*subhi;
if (triclinic == 0) {
sublo = domain->sublo;
subhi = domain->subhi;
} else {
sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda;
}
int m;
double xtmp,ytmp,ztmp;
char *perproc = new char[strlen(file) + 16];
char *ptr = strchr(file,'%');
@ -171,12 +176,15 @@ void ReadRestart::command(int narg, char **arg)
m = 0;
while (m < n) {
xtmp = buf[m+1];
ytmp = buf[m+2];
ztmp = buf[m+3];
if (xtmp >= subxlo && xtmp < subxhi &&
ytmp >= subylo && ytmp < subyhi &&
ztmp >= subzlo && ztmp < subzhi)
x = &buf[m+1];
if (triclinic) {
domain->x2lamda(x,lamda);
coord = lamda;
} else coord = x;
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2])
m += avec->unpack_restart(&buf[m]);
else m += static_cast<int> (buf[m]);
}
@ -520,6 +528,16 @@ void ReadRestart::header()
} else if (flag == 45) {
force->special_coul[3] = read_double();
} else if (flag == 46) {
domain->triclinic = 1;
domain->xy = read_double();
} else if (flag == 47) {
domain->triclinic = 1;
domain->xz = read_double();
} else if (flag == 48) {
domain->triclinic = 1;
domain->yz = read_double();
} else error->all("Invalid flag in header of restart file");
flag = read_int();

18
src/region_block.cpp
View File

@ -28,37 +28,43 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
if (strcmp(arg[2],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
xlo = domain->boxxlo;
if (domain->triclinic == 0) xlo = domain->boxlo[0];
else xlo = domain->boxlo_bound[0];
} else xlo = xscale*atof(arg[2]);
if (strcmp(arg[3],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
xhi = domain->boxxhi;
if (domain->triclinic == 0) xhi = domain->boxhi[0];
else xhi = domain->boxhi_bound[0];
} else xhi = xscale*atof(arg[3]);
if (strcmp(arg[4],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
ylo = domain->boxylo;
if (domain->triclinic == 0) ylo = domain->boxlo[1];
else ylo = domain->boxlo_bound[1];
} else ylo = yscale*atof(arg[4]);
if (strcmp(arg[5],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
yhi = domain->boxyhi;
if (domain->triclinic == 0) yhi = domain->boxhi[1];
else yhi = domain->boxhi_bound[1];
} else yhi = yscale*atof(arg[5]);
if (strcmp(arg[6],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
zlo = domain->boxzlo;
if (domain->triclinic == 0) zlo = domain->boxlo[2];
else zlo = domain->boxlo_bound[2];
} else zlo = zscale*atof(arg[6]);
if (strcmp(arg[7],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
zhi = domain->boxzhi;
if (domain->triclinic == 0) zhi = domain->boxhi[2];
else zhi = domain->boxhi_bound[2];
} else zhi = zscale*atof(arg[7]);
// error check

30
src/region_cylinder.cpp
View File

@ -46,9 +46,18 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[6],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
if (axis == 'x') lo = domain->boxxlo;
if (axis == 'y') lo = domain->boxylo;
if (axis == 'z') lo = domain->boxzlo;
if (axis == 'x') {
if (domain->triclinic == 0) lo = domain->boxlo[0];
else lo = domain->boxlo_bound[0];
}
if (axis == 'y') {
if (domain->triclinic == 0) lo = domain->boxlo[1];
else lo = domain->boxlo_bound[1];
}
if (axis == 'z') {
if (domain->triclinic == 0) lo = domain->boxlo[2];
else lo = domain->boxlo_bound[2];
}
} else {
if (axis == 'x') lo = xscale*atof(arg[6]);
if (axis == 'y') lo = yscale*atof(arg[6]);
@ -58,9 +67,18 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[7],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
if (axis == 'x') hi = domain->boxxhi;
if (axis == 'y') hi = domain->boxyhi;
if (axis == 'z') hi = domain->boxzhi;
if (axis == 'x') {
if (domain->triclinic == 0) hi = domain->boxhi[0];
else hi = domain->boxhi_bound[0];
}
if (axis == 'y') {
if (domain->triclinic == 0) hi = domain->boxhi[1];
else hi = domain->boxhi_bound[1];
}
if (axis == 'z') {
if (domain->triclinic == 0) hi = domain->boxhi[2];
else hi = domain->boxhi_bound[2];
}
} else {
if (axis == 'x') hi = xscale*atof(arg[7]);
if (axis == 'y') hi = yscale*atof(arg[7]);

43
src/region_prism.cpp
View File

@ -36,42 +36,42 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
if (strcmp(arg[2],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
xlo = domain->boxxlo;
xlo = domain->boxlo[0];
} else xlo = xscale*atof(arg[2]);
if (strcmp(arg[3],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
xhi = domain->boxxhi;
xhi = domain->boxhi[0];
} else xhi = xscale*atof(arg[3]);
if (strcmp(arg[4],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
ylo = domain->boxylo;
ylo = domain->boxlo[1];
} else ylo = yscale*atof(arg[4]);
if (strcmp(arg[5],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
yhi = domain->boxyhi;
yhi = domain->boxhi[1];
} else yhi = yscale*atof(arg[5]);
if (strcmp(arg[6],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
zlo = domain->boxzlo;
zlo = domain->boxlo[0];
} else zlo = zscale*atof(arg[6]);
if (strcmp(arg[7],"INF") == 0) {
if (domain->box_exist == 0)
error->all("Cannot use region INF when box does not exist");
zhi = domain->boxzhi;
zhi = domain->boxhi[1];
} else zhi = zscale*atof(arg[7]);
yxshift = xscale*atof(arg[8]);
zxshift = xscale*atof(arg[9]);
zyshift = yscale*atof(arg[10]);
xy = xscale*atof(arg[8]);
xz = xscale*atof(arg[9]);
yz = yscale*atof(arg[10]);
// error check
// prism cannot be 0 thickness in any dim, else inverse blows up
@ -81,13 +81,14 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
// extent of prism
extent_xlo = MIN(xlo,xlo+yxshift);
extent_xlo = MIN(extent_xlo,extent_xlo+zxshift);
extent_xhi = MAX(xhi,xhi+yxshift);
extent_xhi = MAX(extent_xhi,extent_xhi+zxshift);
extent_ylo = MIN(ylo,ylo+zyshift);
extent_yhi = MAX(yhi,yhi+zyshift);
extent_xlo = MIN(xlo,xlo+xy);
extent_xlo = MIN(extent_xlo,extent_xlo+xz);
extent_ylo = MIN(ylo,ylo+yz);
extent_zlo = zlo;
extent_xhi = MAX(xhi,xhi+xy);
extent_xhi = MAX(extent_xhi,extent_xhi+xz);
extent_yhi = MAX(yhi,yhi+yz);
extent_zhi = zhi;
// h = transformation matrix from tilt coords (0-1) to box coords (xyz)
@ -98,27 +99,27 @@ RegPrism::RegPrism(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
// and bottom face of prism is in xy plane
h[0][0] = xhi - xlo;
h[0][1] = yxshift;
h[0][1] = xy;
h[0][2] = xz;
h[1][1] = yhi - ylo;
h[0][2] = zxshift;
h[1][2] = zyshift;
h[1][2] = yz;
h[2][2] = zhi - zlo;
hinv[0][0] = 1.0/h[0][0];
hinv[0][1] = -h[0][1] / (h[0][0]*h[1][1]);
hinv[1][1] = 1.0/h[1][1];
hinv[0][2] = (h[0][1]*h[1][2] - h[0][2]*h[1][1]) / (h[0][0]*h[1][1]*h[2][2]);
hinv[1][1] = 1.0/h[1][1];
hinv[1][2] = -h[1][2] / (h[1][1]*h[2][2]);
hinv[2][2] = 1.0/h[2][2];
}
/* ----------------------------------------------------------------------
check xyz against prism
abc = Hinv * (xyz - xyzlo)
abc = Hinv * (xyz - xyz/lo)
abc = tilt coords (0-1)
Hinv = transformation matrix from box coords to tilt coords
xyz = box coords
xyzlo = lower-left corner of prism
xyz/lo = lower-left corner of prism
------------------------------------------------------------------------- */
int RegPrism::match(double x, double y, double z)

4
src/region_prism.h
View File

@ -19,13 +19,15 @@
namespace LAMMPS_NS {
class RegPrism : public Region {
friend class CreateBox;
public:
RegPrism(class LAMMPS *, int, char **);
int match(double, double, double);
private:
double xlo,xhi,ylo,yhi,zlo,zhi;
double yxshift,zxshift,zyshift;
double xy,xz,yz;
double h[3][3],hinv[3][3];
};

95
src/replicate.cpp
View File

@ -27,7 +27,8 @@
using namespace LAMMPS_NS;
#define LB_FACTOR 1.1
#define MAXATOMS 0x7FFFFFFF
#define MAXATOMS 0x7FFFFFFF
#define EPSILON 1.0e-6
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
@ -92,7 +93,7 @@ void Replicate::command(int narg, char **arg)
// unmap existing atoms via image flags
for (i = 0; i < atom->nlocal; i++)
domain->unmap(atom->x[i][0],atom->x[i][1],atom->x[i][2],atom->image[i]);
domain->unmap(atom->x[i],atom->image[i]);
// communication buffer for all my atom's info
// max_size = largest buffer needed by any proc
@ -157,15 +158,24 @@ void Replicate::command(int narg, char **arg)
// store old simulation box
int triclinic = domain->triclinic;
double old_xprd = domain->xprd;
double old_yprd = domain->yprd;
double old_zprd = domain->zprd;
double old_xy = domain->xy;
double old_xz = domain->xz;
double old_yz = domain->yz;
// setup new simulation box
domain->boxxhi = domain->boxxlo + nx*old_xprd;
domain->boxyhi = domain->boxylo + ny*old_yprd;
domain->boxzhi = domain->boxzlo + nz*old_zprd;
if (triclinic) {
domain->xy *= ny;
domain->xz *= nz;
domain->yz *= nz;
}
// new problem setup using new box boundaries
@ -175,20 +185,12 @@ void Replicate::command(int narg, char **arg)
atom->allocate_type_arrays();
atom->avec->grow(n);
domain->print_box(" ");
domain->set_initial_box();
domain->set_global_box();
comm->set_procs();
domain->set_local_box();
// sub xyz lo/hi = new processor sub-domain
double subxlo = domain->subxlo;
double subxhi = domain->subxhi;
double subylo = domain->subylo;
double subyhi = domain->subyhi;
double subzlo = domain->subzlo;
double subzhi = domain->subzhi;
// copy type arrays to new atom class
if (atom->mass) {
@ -206,12 +208,41 @@ void Replicate::command(int narg, char **arg)
}
}
// set bounds for my proc
// if periodic and I am lo/hi proc, adjust bounds by EPSILON
// insures all replicated atoms will be owned even with round-off
double sublo[3],subhi[3];
if (triclinic == 0) {
sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0];
sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1];
sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2];
} else {
sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0];
sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1];
sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2];
}
if (domain->xperiodic) {
if (comm->myloc[0] == 0) sublo[0] -= EPSILON;
if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += EPSILON;
}
if (domain->yperiodic) {
if (comm->myloc[1] == 0) sublo[1] -= EPSILON;
if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += EPSILON;
}
if (domain->zperiodic) {
if (comm->myloc[2] == 0) sublo[2] -= EPSILON;
if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += EPSILON;
}
// loop over all procs
// if this iteration of loop is me:
// pack my unmapped atom data into buf
// bcast it to all other procs
// for each atom in buf, each proc performs 3d replicate loop:
// xnew,ynew,znew = new replicated position
// performs 3d replicate loop with while loop over atoms in buf
// x = new replicated position
// unpack atom into new atom class from buf if I own it
// adjust tag, mol #, coord, topology info as needed
@ -219,7 +250,8 @@ void Replicate::command(int narg, char **arg)
AtomVec *avec = atom->avec;
int ix,iy,iz,image,atom_offset,mol_offset;
double xnew,ynew,znew;
double x[3],lamda[3];
double *coord;
int tag_enable = atom->tag_enable;
for (int iproc = 0; iproc < nprocs; iproc++) {
@ -234,17 +266,29 @@ void Replicate::command(int narg, char **arg)
for (iy = 0; iy < ny; iy++) {
for (iz = 0; iz < nz; iz++) {
// while loop over one proc's atom list
m = 0;
while (m < n) {
xnew = buf[m+1] + ix*old_xprd;
ynew = buf[m+2] + iy*old_yprd;
znew = buf[m+3] + iz*old_zprd;
image = (512 << 20) | (512 << 10) | 512;
domain->remap(xnew,ynew,znew,image);
if (triclinic == 0) {
x[0] = buf[m+1] + ix*old_xprd;
x[1] = buf[m+2] + iy*old_yprd;
x[2] = buf[m+3] + iz*old_zprd;
} else {
x[0] = buf[m+1] + ix*old_xprd + iy*old_xy + iz*old_xz;
x[1] = buf[m+2] + iy*old_yprd + iz*old_yz;
x[2] = buf[m+3] + iz*old_zprd;
}
domain->remap(x,image);
if (triclinic) {
domain->x2lamda(x,lamda);
coord = lamda;
} else coord = x;
if (xnew >= subxlo && xnew < subxhi &&
ynew >= subylo && ynew < subyhi &&
znew >= subzlo && znew < subzhi) {
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) {
m += avec->unpack_restart(&buf[m]);
@ -254,9 +298,9 @@ void Replicate::command(int narg, char **arg)
else atom_offset = 0;
mol_offset = iz*ny*nx*maxmol + iy*nx*maxmol + ix*maxmol;
atom->x[i][0] = xnew;
atom->x[i][1] = ynew;
atom->x[i][2] = znew;
atom->x[i][0] = x[0];
atom->x[i][1] = x[1];
atom->x[i][2] = x[2];
atom->tag[i] += atom_offset;
atom->image[i] = image;
@ -289,8 +333,7 @@ void Replicate::command(int narg, char **arg)
}
}
} else m += static_cast<int> (buf[m]);
} // end of while loop over one proc's atom list
}
}
}
}

8
src/respa.cpp
View File

@ -321,6 +321,10 @@ void Respa::init()
newton[ilevel] = 1;
}
}
// orthogonal vs triclinic simulation box
triclinic = domain->triclinic;
}
/* ----------------------------------------------------------------------
@ -335,12 +339,14 @@ void Respa::setup()
// acquire ghosts
// build neighbor lists
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
neighbor->build();
neighbor->ncalls = 0;
@ -464,6 +470,7 @@ void Respa::recurse(int ilevel)
int nflag = neighbor->decide();
if (nflag) {
if (modify->n_pre_exchange) modify->pre_exchange();
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
@ -473,6 +480,7 @@ void Respa::recurse(int ilevel)
timer->stamp();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(TIME_COMM);
if (modify->n_pre_neighbor) modify->pre_neighbor();
neighbor->build();

1
src/respa.h
View File

@ -51,6 +51,7 @@ class Respa : public Integrate {
int nfix_virial; // # of fixes that need virial occasionally
int *fix_virial_every; // frequency they require it
int *next_fix_virial; // next timestep they need it
int triclinic; // 0 if domain is orthog, 1 if triclinic
int *newton; // newton flag at each level
class FixRespa *fix_respa; // Fix to store the force level array

2
src/set.cpp
View File

@ -103,11 +103,13 @@ void Set::command(int narg, char **arg)
if (comm->me == 0 && screen) fprintf(screen,"System init for set ...\n");
lmp->init();
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
}
if (comm->me == 0 && screen) fprintf(screen,"Setting atom values ...\n");

28
src/velocity.cpp
View File

@ -283,7 +283,7 @@ void Velocity::create(int narg, char **arg)
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
triple(x[i][0],x[i][1],x[i][2],&vx,&vy,&vz,seed,random);
triple(x[i],&vx,&vy,&vz,seed,random);
if (mass) factor = 1.0/sqrt(mass[type[i]]);
else factor = 1.0/sqrt(rmass[i]);
v[i][0] = vx * factor;
@ -661,31 +661,37 @@ void Velocity::options(int narg, char **arg)
#define IC3 54773
#define IM3 259200
void Velocity::triple(double x, double y, double z,
double *vx, double *vy, double *vz,
void Velocity::triple(double *x, double *vx, double *vy, double *vz,
int seed, RanPark *random)
{
// lamda[3] = fractional position in box
// for triclinic box, just convert to lamda coords
double lamda[3];
if (domain->triclinic == 0) {
lamda[0] = (x[0] - domain->boxlo[0]) / domain->prd[0];
lamda[1] = (x[1] - domain->boxlo[1]) / domain->prd[1];
lamda[2] = (x[2] - domain->boxlo[2]) / domain->prd[2];
} else domain->x2lamda(x,lamda);
// seed 1,2,3 = combination of atom coord in each dim and user-input seed
// map geometric extent into range of each of 3 RNGs
// warm-up each RNG by calling it twice
double fraction;
int seed1,seed2,seed3;
fraction = (x - domain->boxxlo) / domain->xprd;
seed1 = static_cast<int> (fraction * IM1);
seed1 = static_cast<int> (lamda[0] * IM1);
seed1 = (seed1+seed) % IM1;
seed1 = (seed1*IA1+IC1) % IM1;
seed1 = (seed1*IA1+IC1) % IM1;
fraction = (y - domain->boxylo) / domain->yprd;
seed2 = static_cast<int> (fraction * IM2);
seed2 = static_cast<int> (lamda[1] * IM2);
seed2 = (seed2+seed) % IM2;
seed2 = (seed2*IA2+IC2) % IM2;
seed2 = (seed2*IA2+IC2) % IM2;
fraction = (z - domain->boxzlo) / domain->zprd;
seed3 = static_cast<int> (fraction * IM3);
seed3 = static_cast<int> (lamda[2] * IM3);
seed3 = (seed3+seed) % IM3;
seed3 = (seed3*IA3+IC3) % IM3;
seed3 = (seed3*IA3+IC3) % IM3;
@ -693,7 +699,7 @@ void Velocity::triple(double x, double y, double z,
// fraction = 0-1 with giving each dim an equal weighting
// use fraction to reset Park/Miller RNG seed
fraction = 1.0*seed1/(3*IM1) + 1.0*seed2/(3*IM2) + 1.0*seed3/(3*IM3);
double fraction = 1.0*seed1/(3*IM1) + 1.0*seed2/(3*IM2) + 1.0*seed3/(3*IM3);
random->reset(fraction);
// use RNG to set velocities after warming up twice

2
src/velocity.h
View File

@ -41,7 +41,7 @@ class Velocity : protected Pointers {
void zero_rotation();
void options(int, char **);
void triple(double, double, double, double *, double *, double *,
void triple(double *, double *, double *, double *,
int, class RanPark *);
};

8
src/verlet.cpp
View File

@ -103,6 +103,10 @@ void Verlet::init()
pairflag = 1;
if (strcmp(atom->atom_style,"granular") == 0) pairflag = 0;
// orthogonal vs triclinic simulation box
triclinic = domain->triclinic;
// local copies of Update quantities
maxpair = update->maxpair;
@ -121,12 +125,14 @@ void Verlet::setup()
// acquire ghosts
// build neighbor lists
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
neighbor->build();
neighbor->ncalls = 0;
@ -197,6 +203,7 @@ void Verlet::iterate(int n)
timer->stamp(TIME_COMM);
} else {
if (modify->n_pre_exchange) modify->pre_exchange();
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
@ -206,6 +213,7 @@ void Verlet::iterate(int n)
timer->stamp();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(TIME_COMM);
if (modify->n_pre_neighbor) modify->pre_neighbor();
neighbor->build();

1
src/verlet.h
View File

@ -33,6 +33,7 @@ class Verlet : public Integrate {
int nfix_virial; // # of fixes that need virial occasionally
int *fix_virial_every; // frequency they require it
int *next_fix_virial; // next timestep they need it
int triclinic; // 0 if domain is orthog, 1 if triclinic
int pairflag,torqueflag,granflag; // arrays to zero out every step
int maxpair; // local copies of Update quantities

8
src/write_restart.cpp
View File

@ -77,10 +77,12 @@ void WriteRestart::command(int narg, char **arg)
// move atoms to new processors before writing file
// enforce PBC before in case atoms are outside box
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
comm->exchange();
if (domain->triclinic) domain->lamda2x(atom->nlocal);
write(file);
delete [] file;
@ -284,6 +286,12 @@ void WriteRestart::header()
write_double(44,force->special_coul[2]);
write_double(45,force->special_coul[3]);
if (domain->triclinic) {
write_double(46,domain->xy);
write_double(47,domain->xz);
write_double(48,domain->yz);
}
// -1 flag signals end of header
int flag = -1;