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

This commit is contained in:
sjplimp 2010-09-15 22:16:20 +00:00
parent 1d3b37782f
commit ded8b22c19
19 changed files with 909 additions and 747 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

190
src/XTC/dump_xtc.cpp
View File

@ -15,10 +15,11 @@
Contributing authors: Naveen Michaud-Agrawal (Johns Hopkins U)
open-source XDR routines from
Frans van Hoesel (http://md.chem.rug.nl/hoesel)
are also included in this file
Axel Kohlmeyer (UPenn)
are included in this file
Axel Kohlmeyer (Temple U)
port to platforms without XDR support
added support for unwrapped trajectories
support for groups
------------------------------------------------------------------------- */
#include "math.h"
@ -31,6 +32,7 @@
#include "atom.h"
#include "update.h"
#include "group.h"
#include "output.h"
#include "error.h"
#include "memory.h"
@ -52,11 +54,12 @@ int xdr3dfcoord(XDR *, float *, int *, float *);
DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
{
if (narg != 5) error->all("Illegal dump xtc command");
if (igroup != group->find("all")) error->all("Dump xtc must use group all");
if (binary || compressed || multifile || multiproc)
error->all("Invalid dump xtc filename");
size_one = 4;
size_one = 3;
sort_flag = 1;
sortcol = 0;
format_default = NULL;
flush_flag = 0;
unwrap_flag = 0;
@ -64,10 +67,10 @@ DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
// allocate global array for atom coords
natoms = static_cast<int> (atom->natoms);
if (igroup == 0) natoms = static_cast<int> (atom->natoms);
else natoms = static_cast<int> (group->count(igroup));
if (natoms <= 0) error->all("Invalid natoms for dump xtc");
if (atom->tag_consecutive() == 0)
error->all("Atom IDs must be consecutive for dump xtc");
coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords");
// sfactor = conversion of coords to XTC units
@ -77,6 +80,7 @@ DumpXTC::DumpXTC(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
if (strcmp(update->unit_style,"lj") == 0) sfactor = 1.0;
openfile();
nevery_save = 0;
ntotal = 0;
}
@ -94,45 +98,26 @@ DumpXTC::~DumpXTC()
/* ---------------------------------------------------------------------- */
void DumpXTC::init()
void DumpXTC::init_style()
{
if (sort_flag == 0 || sortcol != 0)
error->all("Dump xtc requires sorting by atom ID");
// check that flush_flag is not set since dump::write() will use it
if (flush_flag) error->all("Cannot set dump_modify flush for dump xtc");
}
/* ---------------------------------------------------------------------- */
// check that dump frequency has not changed and is not a variable
int DumpXTC::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"unwrap") == 0) {
if (narg < 2) error->all("Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
else error->all("Illegal dump_modify command");
return 2;
} else if (strcmp(arg[0],"precision") == 0) {
if (narg < 2) error->all("Illegal dump_modify command");
precision = atof(arg[1]);
if ((fabs(precision-10.0) > EPS) && (fabs(precision-100.0) > EPS) &&
(fabs(precision-1000.0) > EPS) && (fabs(precision-10000.0) > EPS) &&
(fabs(precision-100000.0) > EPS) &&
(fabs(precision-1000000.0) > EPS))
error->all("Illegal dump_modify command");
return 2;
}
return 0;
}
int idump;
for (idump = 0; idump < output->ndump; idump++)
if (strcmp(id,output->dump[idump]->id) == 0) break;
if (output->every_dump[idump] == 0)
error->all("Cannot use variable every setting for dump dcd");
/* ----------------------------------------------------------------------
return # of bytes of allocated memory in buf and global coords array
------------------------------------------------------------------------- */
double DumpXTC::memory_usage()
{
double bytes = maxbuf * sizeof(double);
bytes += 3*natoms * sizeof(float);
return bytes;
if (nevery_save == 0) nevery_save = output->every_dump[idump];
else if (nevery_save != output->every_dump[idump])
error->all("Cannot change dump_modify every for dump dcd");
}
/* ---------------------------------------------------------------------- */
@ -151,13 +136,11 @@ void DumpXTC::openfile()
void DumpXTC::write_header(int n)
{
// all procs realloc types & coords if necessary
// all procs realloc coords if total count grew
if (n != natoms) {
memory->sfree(coords);
if (atom->tag_consecutive() == 0)
error->all("Atom IDs must be consecutive for dump xtc");
natoms = n;
memory->sfree(coords);
coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords");
}
@ -201,21 +184,32 @@ void DumpXTC::write_header(int n)
int DumpXTC::count()
{
return atom->nlocal;
if (igroup == 0) return atom->nlocal;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int m = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) m++;
return m;
}
/* ---------------------------------------------------------------------- */
int DumpXTC::pack()
void DumpXTC::pack(int *ids)
{
int m,n;
int *tag = atom->tag;
double **x = atom->x;
int *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// assume group all, so no need to perform mask check
int m = 0;
m = n = 0;
if (unwrap_flag == 1) {
double xprd = domain->xprd;
double yprd = domain->yprd;
@ -224,55 +218,57 @@ int DumpXTC::pack()
double xz = domain->xz;
double yz = domain->yz;
for (int i = 0; i < nlocal; i++) {
int ix = (image[i] & 1023) - 512;
int iy = (image[i] >> 10 & 1023) - 512;
int iz = (image[i] >> 20) - 512;
if (domain->triclinic) {
buf[m++] = tag[i];
buf[m++] = sfactor * (x[i][0] + ix * xprd + iy * xy + iz * xz);
buf[m++] = sfactor * (x[i][1] + iy * yprd + iz * yz);
buf[m++] = sfactor * (x[i][2] + iz * zprd);
} else {
buf[m++] = tag[i];
buf[m++] = sfactor * (x[i][0] + ix * xprd);
buf[m++] = sfactor * (x[i][1] + iy * yprd);
buf[m++] = sfactor * (x[i][2] + iz * zprd);
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
int ix = (image[i] & 1023) - 512;
int iy = (image[i] >> 10 & 1023) - 512;
int iz = (image[i] >> 20) - 512;
if (domain->triclinic) {
buf[m++] = tag[i];
buf[m++] = sfactor * (x[i][0] + ix * xprd + iy * xy + iz * xz);
buf[m++] = sfactor * (x[i][1] + iy * yprd + iz * yz);
buf[m++] = sfactor * (x[i][2] + iz * zprd);
} else {
buf[m++] = tag[i];
buf[m++] = sfactor * (x[i][0] + ix * xprd);
buf[m++] = sfactor * (x[i][1] + iy * yprd);
buf[m++] = sfactor * (x[i][2] + iz * zprd);
}
ids[n++] = tag[i];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
buf[m++] = tag[i];
buf[m++] = sfactor*x[i][0];
buf[m++] = sfactor*x[i][1];
buf[m++] = sfactor*x[i][2];
}
}
return m;
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
buf[m++] = tag[i];
buf[m++] = sfactor*x[i][0];
buf[m++] = sfactor*x[i][1];
buf[m++] = sfactor*x[i][2];
ids[n++] = tag[i];
}
}
}
/* ---------------------------------------------------------------------- */
void DumpXTC::write_data(int n, double *mybuf)
{
float *xyz;
int j,tag;
int j;
// copy buf atom coords into global array
int m = 0;
int k = 3*ntotal;
for (int i = 0; i < n; i++) {
tag = static_cast<int> (mybuf[m]) - 1;
j = 3*tag;
coords[j++] = mybuf[m+1];
coords[j++] = mybuf[m+2];
coords[j] = mybuf[m+3];
m += size_one;
coords[k++] = mybuf[m++];
coords[k++] = mybuf[m++];
coords[k++] = mybuf[m++];
ntotal++;
}
// if last chunk of atoms in this snapshot, write global arrays to file
ntotal += n;
if (ntotal == natoms) {
write_frame();
ntotal = 0;
@ -281,6 +277,40 @@ void DumpXTC::write_data(int n, double *mybuf)
/* ---------------------------------------------------------------------- */
int DumpXTC::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"unwrap") == 0) {
if (narg < 2) error->all("Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
else error->all("Illegal dump_modify command");
return 2;
} else if (strcmp(arg[0],"precision") == 0) {
if (narg < 2) error->all("Illegal dump_modify command");
precision = atof(arg[1]);
if ((fabs(precision-10.0) > EPS) && (fabs(precision-100.0) > EPS) &&
(fabs(precision-1000.0) > EPS) && (fabs(precision-10000.0) > EPS) &&
(fabs(precision-100000.0) > EPS) &&
(fabs(precision-1000000.0) > EPS))
error->all("Illegal dump_modify command");
return 2;
}
return 0;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory in buf and global coords array
------------------------------------------------------------------------- */
double DumpXTC::memory_usage()
{
double bytes = Dump::memory_usage();
bytes += 3*natoms * sizeof(float);
return bytes;
}
/* ---------------------------------------------------------------------- */
void DumpXTC::write_frame()
{
xdr3dfcoord(&xd,coords,&natoms,&precision);

7
src/XTC/dump_xtc.h
View File

@ -35,23 +35,24 @@ class DumpXTC : public Dump {
public:
DumpXTC(class LAMMPS *, int, char**);
~DumpXTC();
void init();
double memory_usage();
private:
int natoms,ntotal;
int nevery_save;
int unwrap_flag; // 1 if atom coords are unwrapped, 0 if no
float precision; // user-adjustable precision setting
float *coords;
double sfactor;
XDR xd;
void init_style();
int modify_param(int, char **);
void openfile();
void write_header(int);
int count();
int pack();
void pack(int *);
void write_data(int, double *);
double memory_usage();
void write_frame();
};

5
src/XTC/xdr_compat.h
View File

@ -60,6 +60,11 @@ extern "C" {
typedef int bool_t;
#if defined(__MINGW32_VERSION)
typedef char * caddr_t;
typedef unsigned int u_int;
#endif
/*
* Aninteger type that is 32 bits wide. Check if int,
* long or short is 32 bits and die if none of them is :-)

383
src/dump.cpp
View File

@ -17,6 +17,7 @@
#include "stdio.h"
#include "dump.h"
#include "atom.h"
#include "irregular.h"
#include "update.h"
#include "domain.h"
#include "group.h"
@ -26,6 +27,17 @@
using namespace LAMMPS_NS;
// allocate space for static class variable
Dump *Dump::dumpptr;
#define BIG 1.0e20
#define IBIG 2147483647
#define EPSILON 1.0e-6
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
/* ---------------------------------------------------------------------- */
Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
@ -56,8 +68,10 @@ Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
sort_flag = 0;
append_flag = 0;
maxbuf = 0;
buf = NULL;
maxbuf = maxsort = maxproc = 0;
buf = bufsort = NULL;
idsort = index = proclist = NULL;
irregular = NULL;
// parse filename for special syntax
// if contains '%', write one file per proc and replace % with proc-ID
@ -108,6 +122,11 @@ Dump::~Dump()
delete [] format_user;
memory->sfree(buf);
memory->sfree(bufsort);
memory->sfree(idsort);
memory->sfree(index);
memory->sfree(proclist);
delete irregular;
// XTC style sets fp to NULL since it closes file in its destructor
@ -124,6 +143,36 @@ Dump::~Dump()
/* ---------------------------------------------------------------------- */
void Dump::init()
{
init_style();
if (!sort_flag) {
memory->sfree(buf);
memory->sfree(bufsort);
memory->sfree(idsort);
memory->sfree(index);
memory->sfree(proclist);
delete irregular;
maxbuf = maxsort = maxproc = 0;
buf = bufsort = NULL;
idsort = index = proclist = NULL;
irregular = NULL;
}
if (sort_flag && sortcol == 0 && atom->tag_enable == 0)
error->all("Cannot use dump sort on atom IDs with no atom IDs defined");
if (sort_flag && sortcol > size_one)
error->all("Dump sort column is invalid");
if (sort_flag && nprocs > 1 && irregular == NULL)
irregular = new Irregular(lmp);
}
/* ---------------------------------------------------------------------- */
void Dump::write()
{
// if file per timestep, open new file
@ -152,11 +201,12 @@ void Dump::write()
}
// nme = # of dump lines this proc will contribute to dump
nme = count();
// ntotal = total # of dump lines
// nmax = max # of dump lines on any proc
int nme = count();
int ntotal,nmax;
if (multiproc) nmax = nme;
else {
@ -169,25 +219,31 @@ void Dump::write()
if (multiproc) write_header(nme);
else write_header(ntotal);
// grow communication buffer if necessary
// pack my data into buf
// if sorting on IDs also request ID list from pack()
// sort buf as needed
if (nmax*size_one > maxbuf) {
maxbuf = nmax*size_one;
if (nmax > maxbuf) {
maxbuf = nmax;
memory->sfree(buf);
buf = (double *) memory->smalloc(maxbuf*sizeof(double),"dump:buf");
buf = (double *)
memory->smalloc(maxbuf*size_one*sizeof(double),"dump:buf");
if (sort_flag && sortcol == 0) {
memory->sfree(ids);
ids = (int *) memory->smalloc(maxbuf*sizeof(int),"dump:ids");
}
}
// pack my data into buf
// me_size = # of quantities in buf
int me_size = pack();
if (sort_flag && sortcol == 0) pack(ids);
else pack(NULL);
if (sort_flag) sort();
// multiproc = 1 = each proc writes own data to own file
// multiproc = 0 = all procs write to one file thru proc 0
// proc 0 pings each proc, receives it's data, writes to file
// all other procs wait for ping, send their data to proc 0
if (multiproc) write_data(me_size/size_one,buf);
if (multiproc) write_data(nme,buf);
else {
int tmp,nlines;
MPI_Status status;
@ -201,7 +257,7 @@ void Dump::write()
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_DOUBLE,&nlines);
nlines /= size_one;
} else nlines = me_size/size_one;
} else nlines = nme;
write_data(nlines,buf);
}
@ -209,7 +265,7 @@ void Dump::write()
} else {
MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
MPI_Rsend(buf,me_size,MPI_DOUBLE,0,0,world);
MPI_Rsend(buf,nme*size_one,MPI_DOUBLE,0,0,world);
}
}
@ -226,6 +282,229 @@ void Dump::write()
}
}
/* ----------------------------------------------------------------------
generic opening of a dump file
ASCII or binary or gzipped
some derived classes override this function
------------------------------------------------------------------------- */
void Dump::openfile()
{
// single file, already opened, so just return
if (singlefile_opened) return;
if (multifile == 0) singlefile_opened = 1;
// if one file per timestep, replace '*' with current timestep
char *filecurrent;
if (multifile == 0) filecurrent = filename;
else {
filecurrent = new char[strlen(filename) + 16];
char *ptr = strchr(filename,'*');
*ptr = '\0';
sprintf(filecurrent,"%s%d%s",filename,update->ntimestep,ptr+1);
*ptr = '*';
}
// open one file on proc 0 or file on every proc
if (me == 0 || multiproc) {
if (compressed) {
#ifdef LAMMPS_GZIP
char gzip[128];
sprintf(gzip,"gzip -6 > %s",filecurrent);
fp = popen(gzip,"w");
#else
error->one("Cannot open gzipped file");
#endif
} else if (binary) {
fp = fopen(filecurrent,"wb");
} else if (append_flag) {
fp = fopen(filecurrent,"a");
} else {
fp = fopen(filecurrent,"w");
}
if (fp == NULL) error->one("Cannot open dump file");
} else fp = NULL;
// delete string with timestep replaced
if (multifile) delete [] filecurrent;
}
/* ----------------------------------------------------------------------
parallel sort of buf across all procs
changes nme, reorders datums in buf, grows buf if necessary
------------------------------------------------------------------------- */
void Dump::sort()
{
int i,iproc;
double value;
// if single proc, swap ptrs to buf,ids <-> bufsort,idsort
if (nprocs == 1) {
if (nme > maxsort) {
maxsort = nme;
memory->sfree(bufsort);
bufsort = (double *)
memory->smalloc(maxsort*size_one*sizeof(double),"dump:bufsort");
memory->sfree(index);
index = (int *) memory->smalloc(maxsort*sizeof(int),"dump:index");
if (sortcol == 0) {
memory->sfree(idsort);
idsort = (int *) memory->smalloc(maxsort*sizeof(int),"dump:idsort");
}
}
double *dptr = buf;
buf = bufsort;
bufsort = dptr;
if (sortcol == 0) {
int *iptr = ids;
ids = idsort;
idsort = iptr;
}
// if multiple procs, exchange datums between procs via irregular
} else {
// grow proclist if necessary
if (nme > maxproc) {
maxproc = nme;
memory->sfree(proclist);
proclist = (int *) memory->smalloc(maxproc*sizeof(int),"dump:proclist");
}
// proclist[i] = which proc Ith datum will be sent to
if (sortcol == 0) {
int min = IBIG;
int max = 0;
for (i = 0; i < nme; i++) {
min = MIN(min,ids[i]);
max = MAX(max,ids[i]);
}
int minall,maxall;
MPI_Allreduce(&min,&minall,1,MPI_INT,MPI_MIN,world);
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
double range = maxall - minall + 0.1;
for (i = 0; i < nme; i++) {
iproc = static_cast<int> ((ids[i]-minall)/range * nprocs);
proclist[i] = iproc;
}
} else {
double min = BIG;
double max = -BIG;
for (i = 0; i < nme; i++) {
value = buf[i*size_one + sortcolm1];
min = MIN(min,value);
max = MAX(max,value);
}
double minall,maxall;
MPI_Allreduce(&min,&minall,1,MPI_DOUBLE,MPI_MIN,world);
MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world);
double range = maxall-minall + EPSILON*(maxall-minall);
if (range == 0.0) range = EPSILON;
for (i = 0; i < nme; i++) {
value = buf[i*size_one + sortcolm1];
iproc = static_cast<int> ((value-minall)/range * nprocs);
proclist[i] = iproc;
}
}
// create comm plan, grow recv bufs if necessary,
// exchange datums, destroy plan
// if sorting on atom IDs, exchange IDs also
nme = irregular->create_data(nme,proclist);
if (nme > maxsort) {
maxsort = nme;
memory->sfree(bufsort);
bufsort = (double *)
memory->smalloc(maxsort*size_one*sizeof(double),"dump:bufsort");
memory->sfree(index);
index = (int *) memory->smalloc(maxsort*sizeof(int),"dump:index");
if (sortcol == 0) {
memory->sfree(idsort);
idsort = (int *) memory->smalloc(maxsort*sizeof(int),"dump:idsort");
}
}
irregular->exchange_data((char *) buf,size_one*sizeof(double),
(char *) bufsort);
if (sortcol == 0)
irregular->exchange_data((char *) ids,sizeof(int),(char *) idsort);
irregular->destroy_data();
}
// quicksort indices using IDs or buf column as comparator
dumpptr = this;
for (i = 0; i < nme; i++) index[i] = i;
if (sortcol == 0) qsort(index,nme,sizeof(int),idcompare);
else qsort(index,nme,sizeof(int),bufcompare);
// copy data from bufsort to buf using index
if (nme > maxbuf) {
maxbuf = nme;
memory->sfree(buf);
buf = (double *)
memory->smalloc(maxbuf*size_one*sizeof(double),"dump:buf");
}
int nbytes = size_one*sizeof(double);
for (i = 0; i < nme; i++)
memcpy(&buf[i*size_one],&bufsort[index[i]*size_one],nbytes);
}
/* ----------------------------------------------------------------------
compare two atom IDs in sort data structure
called via qsort_r in sort() method
is a static method so access sort data structure via ptr
------------------------------------------------------------------------- */
int Dump::idcompare(const void *pi, const void *pj)
{
int *idsort = dumpptr->idsort;
int i = *((int *) pi);
int j = *((int *) pj);
if (idsort[i] < idsort[j]) return -1;
if (idsort[i] > idsort[j]) return 1;
return 0;
}
/* ----------------------------------------------------------------------
compare two buffer quantities in sort data structure with size_one stride
called via qsort_r in sort() method
is a static method so access sort data structure via ptr
------------------------------------------------------------------------- */
int Dump::bufcompare(const void *pi, const void *pj)
{
double *bufsort = dumpptr->bufsort;
int size_one = dumpptr->size_one;
int sortcolm1 = dumpptr->sortcolm1;
int i = *((int *) pi)*size_one + sortcolm1;
int j = *((int *) pj)*size_one + sortcolm1;
if (bufsort[i] < bufsort[j]) return -1;
if (bufsort[i] > bufsort[j]) return 1;
return 0;
}
/* ----------------------------------------------------------------------
process params common to all dumps here
if unknown param, call modify_param specific to the dump
@ -285,9 +564,13 @@ void Dump::modify_params(int narg, char **arg)
iarg += 2;
} else if (strcmp(arg[iarg],"sort") == 0) {
if (iarg+2 > narg) error->all("Illegal dump_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) sort_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) sort_flag = 0;
else error->all("Illegal dump_modify command");
if (strcmp(arg[iarg+1],"off") == 0) sort_flag = 0;
else {
sort_flag = 1;
sortcol = atoi(arg[iarg+1]);
if (sortcol < 0) error->all("Illegal dump_modify command");
sortcolm1 = sortcol - 1;
}
iarg += 2;
} else {
int n = modify_param(narg-iarg,&arg[iarg]);
@ -298,63 +581,19 @@ void Dump::modify_params(int narg, char **arg)
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory in buf
return # of bytes of allocated memory
------------------------------------------------------------------------- */
double Dump::memory_usage()
{
double bytes = maxbuf * sizeof(double);
double bytes = maxbuf*size_one * sizeof(double);
if (sort_flag) {
if (sortcol == 0) bytes += maxbuf * sizeof(int); // ids
bytes += maxsort*size_one * sizeof(double); // bufsort
bytes += maxsort * sizeof(int); // index
if (sortcol == 0) bytes += maxsort * sizeof(int); // idsort
bytes += maxproc * sizeof(int); // proclist
if (irregular) bytes += irregular->memory_usage();
}
return bytes;
}
/* ----------------------------------------------------------------------
generic opening of a dump file
ASCII or binary or gzipped
some derived classes override this function
------------------------------------------------------------------------- */
void Dump::openfile()
{
// single file, already opened, so just return
if (singlefile_opened) return;
if (multifile == 0) singlefile_opened = 1;
// if one file per timestep, replace '*' with current timestep
char *filecurrent;
if (multifile == 0) filecurrent = filename;
else {
filecurrent = new char[strlen(filename) + 16];
char *ptr = strchr(filename,'*');
*ptr = '\0';
sprintf(filecurrent,"%s%d%s",filename,update->ntimestep,ptr+1);
*ptr = '*';
}
// open one file on proc 0 or file on every proc
if (me == 0 || multiproc) {
if (compressed) {
#ifdef LAMMPS_GZIP
char gzip[128];
sprintf(gzip,"gzip -6 > %s",filecurrent);
fp = popen(gzip,"w");
#else
error->one("Cannot open gzipped file");
#endif
} else if (binary) {
fp = fopen(filecurrent,"wb");
} else if (append_flag) {
fp = fopen(filecurrent,"a");
} else {
fp = fopen(filecurrent,"w");
}
if (fp == NULL) error->one("Cannot open dump file");
} else fp = NULL;
// delete string with timestep replaced
if (multifile) delete [] filecurrent;
}

51
src/dump.h
View File

@ -24,6 +24,22 @@ class Dump : protected Pointers {
char *id; // user-defined name of Dump
char *style; // style of Dump
int igroup,groupbit; // group that Dump is performed on
int first_flag; // 0 if no initial dump, 1 if yes initial dump
int clearstep; // 1 if dump invokes computes, 0 if not
// static variable across all Dump objects
static Dump *dumpptr; // holds a ptr to Dump currently being used
Dump(class LAMMPS *, int, char **);
virtual ~Dump();
void init();
void write();
void modify_params(int, char **);
virtual double memory_usage();
protected:
int me,nprocs; // proc info
char *filename; // user-specified file
@ -33,41 +49,46 @@ class Dump : protected Pointers {
int multiproc; // 0 = proc 0 writes for all, 1 = one file/proc
int header_flag; // 0 = item, 2 = xyz
int first_flag; // 0 if no initial dump, 1 if yes initial dump
int flush_flag; // 0 if no flush, 1 if flush every dump
int sort_flag; // 1 if write in sorted order, 0 if not
int sort_flag; // 1 if sorted output
int append_flag; // 1 if open file in append mode, 0 if not
int singlefile_opened; // 1 = one big file, already opened, else 0
int clearstep; // 1 if dump invokes computes, 0 if not
int sortcol; // 0 to sort on ID, 1-N on columns
int sortcolm1; // sortcol - 1
char *format_default; // default format string
char *format_user; // format string set by user
char *format; // format string for the file write
double *buf; // memory for atom quantities
int maxbuf; // size of buf
FILE *fp; // file to write dump to
int size_one; // # of quantities for one atom
int nme; // # of atoms in this dump from me
Dump(class LAMMPS *, int, char **);
virtual ~Dump();
virtual void init() {}
void write();
void modify_params(int, char **);
virtual double memory_usage();
protected:
double boxxlo,boxxhi; // local copies of domain values
double boxylo,boxyhi; // lo/hi are bounding box for triclinic
double boxzlo,boxzhi;
double boxxy,boxxz,boxyz;
int maxbuf; // size of buf and ids
int maxsort; // size of bufsort, idsort, index
int maxproc; // size of proclist
double *buf; // memory for atom quantities
int *ids; // list of atom IDs, if sorting on IDs
double *bufsort;
int *idsort,*index,*proclist;
class Irregular *irregular;
virtual void init_style() = 0;
virtual void openfile();
virtual int modify_param(int, char **) {return 0;}
virtual void write_header(int) = 0;
virtual int count() = 0;
virtual int pack() = 0;
virtual void pack(int *) = 0;
virtual void write_data(int, double *) = 0;
void sort();
static int idcompare(const void *, const void *);
static int bufcompare(const void *, const void *);
};
}

60
src/dump_atom.cpp
View File

@ -34,7 +34,7 @@ DumpAtom::DumpAtom(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
/* ---------------------------------------------------------------------- */
void DumpAtom::init()
void DumpAtom::init_style()
{
if (image_flag == 0) size_one = 5;
else size_one = 8;
@ -146,9 +146,9 @@ int DumpAtom::count()
/* ---------------------------------------------------------------------- */
int DumpAtom::pack()
void DumpAtom::pack(int *ids)
{
return (this->*pack_choice)();
(this->*pack_choice)(ids);
}
/* ---------------------------------------------------------------------- */
@ -233,8 +233,10 @@ void DumpAtom::header_item_triclinic(int ndump)
/* ---------------------------------------------------------------------- */
int DumpAtom::pack_scale_image()
void DumpAtom::pack_scale_image(int *ids)
{
int m,n;
int *tag = atom->tag;
int *type = atom->type;
int *image = atom->image;
@ -246,7 +248,7 @@ int DumpAtom::pack_scale_image()
double invyprd = 1.0/domain->yprd;
double invzprd = 1.0/domain->zprd;
int m = 0;
m = n = 00;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
buf[m++] = tag[i];
@ -257,15 +259,16 @@ int DumpAtom::pack_scale_image()
buf[m++] = (image[i] & 1023) - 512;
buf[m++] = (image[i] >> 10 & 1023) - 512;
buf[m++] = (image[i] >> 20) - 512;
if (ids) ids[n++] = tag[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
int DumpAtom::pack_scale_noimage()
void DumpAtom::pack_scale_noimage(int *ids)
{
int m,n;
int *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
@ -276,7 +279,7 @@ int DumpAtom::pack_scale_noimage()
double invyprd = 1.0/domain->yprd;
double invzprd = 1.0/domain->zprd;
int m = 0;
m = n = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
buf[m++] = tag[i];
@ -284,15 +287,16 @@ int DumpAtom::pack_scale_noimage()
buf[m++] = (x[i][0] - boxxlo) * invxprd;
buf[m++] = (x[i][1] - boxylo) * invyprd;
buf[m++] = (x[i][2] - boxzlo) * invzprd;
if (ids) ids[n++] = tag[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
int DumpAtom::pack_scale_image_triclinic()
void DumpAtom::pack_scale_image_triclinic(int *ids)
{
int m,n;
int *tag = atom->tag;
int *type = atom->type;
int *image = atom->image;
@ -302,7 +306,7 @@ int DumpAtom::pack_scale_image_triclinic()
double lamda[3];
int m = 0;
m = n = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
buf[m++] = tag[i];
@ -314,15 +318,16 @@ int DumpAtom::pack_scale_image_triclinic()
buf[m++] = (image[i] & 1023) - 512;
buf[m++] = (image[i] >> 10 & 1023) - 512;
buf[m++] = (image[i] >> 20) - 512;
if (ids) ids[n++] = tag[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
int DumpAtom::pack_scale_noimage_triclinic()
void DumpAtom::pack_scale_noimage_triclinic(int *ids)
{
int m,n;
int *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
@ -331,7 +336,7 @@ int DumpAtom::pack_scale_noimage_triclinic()
double lamda[3];
int m = 0;
m = n = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
buf[m++] = tag[i];
@ -340,15 +345,16 @@ int DumpAtom::pack_scale_noimage_triclinic()
buf[m++] = lamda[0];
buf[m++] = lamda[1];
buf[m++] = lamda[2];
if (ids) ids[n++] = tag[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
int DumpAtom::pack_noscale_image()
void DumpAtom::pack_noscale_image(int *ids)
{
int m,n;
int *tag = atom->tag;
int *type = atom->type;
int *image = atom->image;
@ -356,7 +362,7 @@ int DumpAtom::pack_noscale_image()
double **x = atom->x;
int nlocal = atom->nlocal;
int m = 0;
m = n = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
buf[m++] = tag[i];
@ -367,22 +373,23 @@ int DumpAtom::pack_noscale_image()
buf[m++] = (image[i] & 1023) - 512;
buf[m++] = (image[i] >> 10 & 1023) - 512;
buf[m++] = (image[i] >> 20) - 512;
if (ids) ids[n++] = tag[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
int DumpAtom::pack_noscale_noimage()
void DumpAtom::pack_noscale_noimage(int *ids)
{
int m,n;
int *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
double **x = atom->x;
int nlocal = atom->nlocal;
int m = 0;
m = n = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
buf[m++] = tag[i];
@ -390,9 +397,8 @@ int DumpAtom::pack_noscale_noimage()
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
if (ids) ids[n++] = tag[i];
}
return m;
}
/* ---------------------------------------------------------------------- */

18
src/dump_atom.h
View File

@ -27,7 +27,6 @@ namespace LAMMPS_NS {
class DumpAtom : public Dump {
public:
DumpAtom(LAMMPS *, int, char**);
void init();
private:
int scale_flag; // 1 if atom coords are scaled, 0 if no
@ -35,10 +34,11 @@ class DumpAtom : public Dump {
char *columns; // column labels
void init_style();
int modify_param(int, char **);
void write_header(int);
int count();
int pack();
void pack(int *);
void write_data(int, double *);
typedef void (DumpAtom::*FnPtrHeader)(int);
@ -48,14 +48,14 @@ class DumpAtom : public Dump {
void header_item(int);
void header_item_triclinic(int);
typedef int (DumpAtom::*FnPtrPack)();
typedef void (DumpAtom::*FnPtrPack)(int *);
FnPtrPack pack_choice; // ptr to pack functions
int pack_scale_image();
int pack_scale_noimage();
int pack_noscale_image();
int pack_noscale_noimage();
int pack_scale_image_triclinic();
int pack_scale_noimage_triclinic();
void pack_scale_image(int *);
void pack_scale_noimage(int *);
void pack_noscale_image(int *);
void pack_noscale_noimage(int *);
void pack_scale_image_triclinic(int *);
void pack_scale_noimage_triclinic(int *);
typedef void (DumpAtom::*FnPtrData)(int, double *);
FnPtrData write_choice; // ptr to write data functions

166
src/dump_cfg.cpp
View File

@ -49,9 +49,6 @@ DumpCFG::DumpCFG(LAMMPS *lmp, int narg, char **arg) :
// arrays for data rearrangement
recvcounts = new int[nprocs];
displs = new int[nprocs];
tags = NULL;
rbuf = NULL;
nchosen = nlines = 0;
@ -103,10 +100,6 @@ DumpCFG::~DumpCFG()
delete [] typenames;
}
delete [] recvcounts;
delete [] displs;
if (tags) memory->sfree(tags);
if (rbuf) memory->destroy_2d_double_array(rbuf);
if (auxname) {
@ -117,7 +110,7 @@ DumpCFG::~DumpCFG()
/* ---------------------------------------------------------------------- */
void DumpCFG::init()
void DumpCFG::init_style()
{
if (multifile == 0)
error->all("Dump in CFG format requires one snapshot per file");
@ -222,10 +215,6 @@ void DumpCFG::write_header(int n)
if (rbuf) memory->destroy_2d_double_array(rbuf);
rbuf = memory->create_2d_double_array(nchosen,size_one,"dump:rbuf");
}
// create a sorted list of atom IDs on writing proc(s) if necessary
if (sort_flag) create_sorted_tags();
}
/* ----------------------------------------------------------------------
@ -239,147 +228,40 @@ void DumpCFG::write_data(int n, double *mybuf)
int tag_i,index;
double *mass = atom->mass;
// if sort flag is off, transfer data in buf to rbuf directly
// transfer data from buf to rbuf
// if write by proc 0, transfer chunk by chunk
if (sort_flag == 0) {
for (i = 0, m = 0; i < n; i++) {
for (j = 0; j < size_one; j++)
rbuf[nlines][j] = mybuf[m++];
nlines++;
}
for (i = 0, m = 0; i < n; i++) {
for (j = 0; j < size_one; j++)
rbuf[nlines][j] = mybuf[m++];
nlines++;
}
// write data lines in rbuf to file after transfer is done
if (nlines == nchosen) {
for (itype = 1; itype <= ntypes; itype++) {
for (i = 0; i < nchosen; i++)
if (rbuf[i][1] == itype) break;
if (i < nchosen) {
fprintf(fp,"%g\n",mass[itype]);
fprintf(fp,"%s\n",typenames[itype]);
for (; i < nchosen; i++) {
if (rbuf[i][1] == itype) {
for (j = 2; j < size_one; j++) {
if (vtype[j] == INT)
fprintf(fp,vformat[j],static_cast<int> (rbuf[i][j]));
else fprintf(fp,vformat[j],rbuf[i][j]);
}
fprintf(fp,"\n");
// write data lines in rbuf to file after transfer is done
if (nlines == nchosen) {
for (itype = 1; itype <= ntypes; itype++) {
for (i = 0; i < nchosen; i++)
if (rbuf[i][1] == itype) break;
if (i < nchosen) {
fprintf(fp,"%g\n",mass[itype]);
fprintf(fp,"%s\n",typenames[itype]);
for (; i < nchosen; i++) {
if (rbuf[i][1] == itype) {
for (j = 2; j < size_one; j++) {
if (vtype[j] == INT)
fprintf(fp,vformat[j],static_cast<int> (rbuf[i][j]));
else fprintf(fp,vformat[j],rbuf[i][j]);
}
fprintf(fp,"\n");
}
}
}
nlines = 0;
}
// if sort flag is on, transfer data lines in buf to rbuf &
// rearrange them in a sorted order
// if write by proc 0, transfer chunk by chunk
} else {
// sort data lines:
// find the index of the atom ID of each data line within the
// sorted atom IDs array
// transfer the data line to rbuf with its location according
// to this index
for (i = 0, m = 0; i < n; i++) {
tag_i = static_cast<int>(mybuf[m]);
for (j = 0; j < nchosen; j++) {
if (tag_i == tags[j]) {
index = j;
break;
}
}
for (j = 0; j < size_one; j++)
rbuf[index][j] = mybuf[m++];
}
// write data lines in rbuf to file after transfer is done
nlines += n;
if (nlines == nchosen) {
for (itype = 1; itype <= ntypes; itype++) {
for (i = 0; i < nchosen; i++)
if (rbuf[i][1] == itype) break;
if (i < nchosen) {
fprintf(fp,"%g\n",mass[itype]);
fprintf(fp,"%s\n",typenames[itype]);
for (; i < nchosen; i++) {
if (rbuf[i][1] == itype) {
for (j = 2; j < size_one; j++) {
if (vtype[j] == INT)
fprintf(fp,vformat[j],static_cast<int>(rbuf[i][j]));
else fprintf(fp,vformat[j],rbuf[i][j]);
}
fprintf(fp,"\n");
}
}
}
}
nlines = 0;
}
nlines = 0;
}
}
/* ----------------------------------------------------------------------
create a sorted list of atom IDs on writing proc(s)
------------------------------------------------------------------------- */
void DumpCFG::create_sorted_tags()
{
int index;
int *mytags;
int *tag = atom->tag;
int nlocal = atom->nlocal;
index = 0;
// if write by multiproc, each proc has its own atom IDs collected
if (multiproc) {
if (tags) memory->sfree(tags);
tags = (int *) memory->smalloc(nchosen*sizeof(int),"dump:tags");
for (int i = 0; i < nlocal; i++)
if (choose[i]) tags[index++] = tag[i];
// if write by proc 0, gather atom IDs of all the chosen atoms on proc 0
} else {
mytags = (int *) memory->smalloc(nmine*sizeof(int),"dump:mytags");
for (int i = 0; i < nlocal; i++)
if (choose[i]) mytags[index++] = tag[i];
MPI_Gather(&nmine,1,MPI_INT,recvcounts,1,MPI_INT,0,world);
if (me == 0) {
if (tags) memory->sfree(tags);
tags = (int *) memory->smalloc(nchosen*sizeof(int),"dump:tags");
displs[0] = 0;
for (int iproc = 1; iproc < nprocs; iproc++)
displs[iproc] = displs[iproc-1] + recvcounts[iproc-1];
}
MPI_Gatherv(mytags,nmine,MPI_INT,tags,recvcounts,displs,MPI_INT,0,world);
memory->sfree(mytags);
}
if (multiproc || me == 0) qsort(tags,nchosen,sizeof(int),tag_compare);
}
/* ----------------------------------------------------------------------
compare function for qsort
------------------------------------------------------------------------- */
int DumpCFG::tag_compare(const void *itag, const void *jtag)
{
int tag_i = *((int *) itag);
int tag_j = *((int *) jtag);
if (tag_i < tag_j) return -1;
else if (tag_i > tag_j) return 1;
return 0;
}
/* ---------------------------------------------------------------------- */
int DumpCFG::modify_param2(int narg, char **arg)

14
src/dump_cfg.h
View File

@ -28,10 +28,6 @@ class DumpCFG : public DumpCustom {
public:
DumpCFG(class LAMMPS *, int, char **);
~DumpCFG();
void init();
void write_header(int);
void write_data(int, double *);
int modify_param2(int, char **);
private:
int ntypes; // # of atom types
@ -39,13 +35,13 @@ class DumpCFG : public DumpCustom {
char **auxname; // name strings of auxiliary properties
int nchosen; // # of lines to be written on a writing proc
int nlines; // # of lines transferred from buf to rbuf
int *recvcounts; // # of data (tag info) received from all procs
int *displs; // displacement of data (tag info) within buffer
int *tags; // atom IDs of all the chosen atoms collected
double **rbuf; // buf of data lines for data lines rearrangement
void create_sorted_tags();
static int tag_compare(const void *, const void *);
void init_style();
void write_header(int);
void write_data(int, double *);
int modify_param2(int, char **);
};
}

16
src/dump_custom.cpp
View File

@ -159,7 +159,7 @@ DumpCustom::~DumpCustom()
/* ---------------------------------------------------------------------- */
void DumpCustom::init()
void DumpCustom::init_style()
{
delete [] format;
char *str;
@ -779,10 +779,16 @@ int DumpCustom::count()
/* ---------------------------------------------------------------------- */
int DumpCustom::pack()
void DumpCustom::pack(int *ids)
{
for (int n = 0; n < size_one; n++) (this->*pack_choice[n])(n);
return nmine*size_one;
if (ids) {
int *tag = atom->tag;
int nlocal = atom->nlocal;
int n = 0;
for (int i = 0; i < nlocal; i++)
if (choose[i]) ids[n++] = tag[i];
}
}
/* ---------------------------------------------------------------------- */
@ -1437,7 +1443,7 @@ int DumpCustom::modify_param(int narg, char **arg)
nthresh++;
return 4;
// pass along params to child class
// pass along unknown params to child class
} else {
int n = modify_param2(narg,arg);
@ -1453,7 +1459,7 @@ int DumpCustom::modify_param(int narg, char **arg)
double DumpCustom::memory_usage()
{
double bytes = maxbuf * sizeof(double);
double bytes = Dump::memory_usage();
bytes += maxlocal * sizeof(int);
bytes += maxlocal * sizeof(double);
bytes += maxlocal * nvariable * sizeof(double);

6
src/dump_custom.h
View File

@ -28,8 +28,6 @@ class DumpCustom : public Dump {
public:
DumpCustom(class LAMMPS *, int, char **);
virtual ~DumpCustom();
virtual void init();
double memory_usage();
protected:
int nevery; // dump frequency to check Fix against
@ -71,10 +69,12 @@ class DumpCustom : public Dump {
// private methods
virtual void init_style();
virtual void write_header(int);
int count();
int pack();
void pack(int *);
virtual void write_data(int, double *);
double memory_usage();
void parse_fields(int, char **);
int add_compute(char *);

145
src/dump_dcd.cpp
View File

@ -13,6 +13,7 @@
/* ----------------------------------------------------------------------
Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
Axel Kohlmeyer (Temple U), support for groups
------------------------------------------------------------------------- */
#include "math.h"
@ -54,20 +55,22 @@ static inline void fwrite_int32(FILE* fd, uint32_t i)
DumpDCD::DumpDCD(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
{
if (narg != 5) error->all("Illegal dump dcd command");
if (igroup != group->find("all")) error->all("Dump dcd must use group all");
if (binary || compressed || multifile || multiproc)
error->all("Invalid dump dcd filename");
size_one = 4;
size_one = 3;
sort_flag = 1;
sortcol = 0;
unwrap_flag = 0;
format_default = NULL;
// allocate global array for atom coords
natoms = static_cast<int> (atom->natoms);
if (igroup == 0) natoms = static_cast<int> (atom->natoms);
else natoms = static_cast<int> (group->count(igroup));
if (natoms <= 0) error->all("Invalid natoms for dump dcd");
if (atom->tag_consecutive() == 0)
error->all("Atom IDs must be consecutive for dump dcd");
coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords");
xf = &coords[0*natoms];
yf = &coords[1*natoms];
@ -88,8 +91,11 @@ DumpDCD::~DumpDCD()
/* ---------------------------------------------------------------------- */
void DumpDCD::init()
void DumpDCD::init_style()
{
if (sort_flag == 0 || sortcol != 0)
error->all("Dump dcd requires sorting by atom ID");
// check that dump frequency has not changed and is not a variable
int idump;
@ -105,31 +111,6 @@ void DumpDCD::init()
/* ---------------------------------------------------------------------- */
int DumpDCD::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"unwrap") == 0) {
if (narg < 2) error->all("Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
else error->all("Illegal dump_modify command");
return 2;
}
return 0;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory in buf and global coords array
------------------------------------------------------------------------- */
double DumpDCD::memory_usage()
{
double bytes = maxbuf * sizeof(double);
bytes += 3*natoms * sizeof(float);
return bytes;
}
/* ---------------------------------------------------------------------- */
void DumpDCD::openfile()
{
if (me == 0) {
@ -193,21 +174,30 @@ void DumpDCD::write_header(int n)
int DumpDCD::count()
{
return atom->nlocal;
if (igroup == 0) return atom->nlocal;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int m = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) m++;
return m;
}
/* ---------------------------------------------------------------------- */
int DumpDCD::pack()
void DumpDCD::pack(int *ids)
{
int m,n;
int *tag = atom->tag;
double **x = atom->x;
int *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// assume group all, so no need to perform mask check
int m = 0;
m = n = 0;
if (unwrap_flag) {
double xprd = domain->xprd;
double yprd = domain->yprd;
@ -217,53 +207,51 @@ int DumpDCD::pack()
double yz = domain->yz;
for (int i = 0; i < nlocal; i++) {
int ix = (image[i] & 1023) - 512;
int iy = (image[i] >> 10 & 1023) - 512;
int iz = (image[i] >> 20) - 512;
if (mask[i] & groupbit) {
int ix = (image[i] & 1023) - 512;
int iy = (image[i] >> 10 & 1023) - 512;
int iz = (image[i] >> 20) - 512;
if (domain->triclinic) {
buf[m++] = tag[i];
buf[m++] = x[i][0] + ix * xprd + iy * xy + iz * xz;
buf[m++] = x[i][1] + iy * yprd + iz * yz;
buf[m++] = x[i][2] + iz * zprd;
} else {
buf[m++] = tag[i];
buf[m++] = x[i][0] + ix * xprd;
buf[m++] = x[i][1] + iy * yprd;
buf[m++] = x[i][2] + iz * zprd;
if (domain->triclinic) {
buf[m++] = x[i][0] + ix * xprd + iy * xy + iz * xz;
buf[m++] = x[i][1] + iy * yprd + iz * yz;
buf[m++] = x[i][2] + iz * zprd;
} else {
buf[m++] = x[i][0] + ix * xprd;
buf[m++] = x[i][1] + iy * yprd;
buf[m++] = x[i][2] + iz * zprd;
}
ids[n++] = tag[i];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
buf[m++] = tag[i];
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
}
}
return m;
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
ids[n++] = tag[i];
}
}
}
/* ---------------------------------------------------------------------- */
void DumpDCD::write_data(int n, double *mybuf)
{
// spread buf atom coords into global arrays
// copy buf atom coords into 3 global arrays
int tag;
int m = 0;
for (int i = 0; i < n; i++) {
tag = static_cast<int> (mybuf[m]) - 1;
xf[tag] = mybuf[m+1];
yf[tag] = mybuf[m+2];
zf[tag] = mybuf[m+3];
m += size_one;
xf[ntotal] = mybuf[m++];
yf[ntotal] = mybuf[m++];
zf[ntotal] = mybuf[m++];
ntotal++;
}
// if last chunk of atoms in this snapshot, write global arrays to file
ntotal += n;
if (ntotal == natoms) {
write_frame();
ntotal = 0;
@ -272,6 +260,31 @@ void DumpDCD::write_data(int n, double *mybuf)
/* ---------------------------------------------------------------------- */
int DumpDCD::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"unwrap") == 0) {
if (narg < 2) error->all("Illegal dump_modify command");
if (strcmp(arg[1],"yes") == 0) unwrap_flag = 1;
else if (strcmp(arg[1],"no") == 0) unwrap_flag = 0;
else error->all("Illegal dump_modify command");
return 2;
}
return 0;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory in buf and global coords array
------------------------------------------------------------------------- */
double DumpDCD::memory_usage()
{
double bytes = Dump::memory_usage();
bytes += 3*natoms * sizeof(float);
return bytes;
}
/* ---------------------------------------------------------------------- */
void DumpDCD::write_frame()
{
// write coords

8
src/dump_dcd.h
View File

@ -30,21 +30,21 @@ class DumpDCD : public Dump {
public:
DumpDCD(LAMMPS *, int, char**);
~DumpDCD();
void init();
double memory_usage();
private:
int natoms,ntotal,headerflag,nevery_save,nframes;
float *coords,*xf,*yf,*zf;
int unwrap_flag; // 1 if atom coords are unwrapped, 0 if no
void init_style();
void openfile();
void write_header(int);
int count();
int pack();
void pack(int *);
void write_data(int, double *);
int modify_param(int, char **);
double memory_usage();
void write_frame();
void write_dcd_header(const char *);
};

8
src/dump_local.cpp
View File

@ -120,8 +120,11 @@ DumpLocal::~DumpLocal()
/* ---------------------------------------------------------------------- */
void DumpLocal::init()
void DumpLocal::init_style()
{
if (sort_flag && sortcol == 0)
error->all("Dump local cannot sort by atom ID");
delete [] format;
char *str;
if (format_user) str = format_user;
@ -237,10 +240,9 @@ int DumpLocal::count()
/* ---------------------------------------------------------------------- */
int DumpLocal::pack()
void DumpLocal::pack(int *dummy)
{
for (int n = 0; n < size_one; n++) (this->*pack_choice[n])(n);
return nmine*size_one;
}
/* ---------------------------------------------------------------------- */

4
src/dump_local.h
View File

@ -28,7 +28,6 @@ class DumpLocal : public Dump {
public:
DumpLocal(LAMMPS *, int, char **);
~DumpLocal();
void init();
private:
int nevery; // dump frequency to check Fix against
@ -54,10 +53,11 @@ class DumpLocal : public Dump {
char **id_fix; // their IDs
class Fix **fix; // list of ptrs to the Fix objects
void init_style();
int modify_param(int, char **);
void write_header(int);
int count();
int pack();
void pack(int *);
void write_data(int, double *);
void parse_fields(int, char **);

113
src/dump_xyz.cpp
View File

@ -28,39 +28,18 @@ DumpXYZ::DumpXYZ(LAMMPS *lmp, int narg, char **arg) : Dump(lmp, narg, arg)
if (binary || multiproc) error->all("Invalid dump xyz filename");
size_one = 5;
sort_flag = 1;
sortcol = 0;
char *str = (char *) "%d %g %g %g";
int n = strlen(str) + 1;
format_default = new char[n];
strcpy(format_default,str);
// allocate global array for atom coords if group is all
if (igroup == 0) {
natoms = static_cast<int> (atom->natoms);
if (natoms <= 0) error->all("Invalid natoms for dump xyz");
if (atom->tag_consecutive() == 0)
error->all("Atom IDs must be consecutive for dump xyz");
types = (int *) memory->smalloc(natoms*sizeof(int),"dump:types");
coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords");
}
ntotal = 0;
}
/* ---------------------------------------------------------------------- */
DumpXYZ::~DumpXYZ()
{
if (igroup == 0) {
memory->sfree(types);
memory->sfree(coords);
}
}
/* ---------------------------------------------------------------------- */
void DumpXYZ::init()
void DumpXYZ::init_style()
{
delete [] format;
char *str;
@ -77,41 +56,10 @@ void DumpXYZ::init()
if (multifile == 0) openfile();
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory in buf and global coords array
------------------------------------------------------------------------- */
double DumpXYZ::memory_usage()
{
double bytes = maxbuf * sizeof(double);
if (igroup == 0) {
bytes += natoms * sizeof(int);
bytes += 3*natoms * sizeof(float);
}
return bytes;
}
/* ---------------------------------------------------------------------- */
void DumpXYZ::write_header(int n)
{
// check bounds of atom IDs if group is all
// all procs realloc types & coords if necessary
if (igroup == 0) {
if (atom->tag_consecutive() == 0)
error->all("Atom IDs must be consecutive for dump xyz all");
if (n != natoms) {
memory->sfree(types);
memory->sfree(coords);
natoms = n;
types = (int *) memory->smalloc(natoms*sizeof(int),"dump:types");
coords = (float *) memory->smalloc(3*natoms*sizeof(float),"dump:coords");
}
}
// only proc 0 writes header
if (me == 0) {
fprintf(fp,"%d\n",n);
fprintf(fp,"Atoms\n");
@ -135,15 +83,17 @@ int DumpXYZ::count()
/* ---------------------------------------------------------------------- */
int DumpXYZ::pack()
void DumpXYZ::pack(int *ids)
{
int m,n;
int *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
double **x = atom->x;
int nlocal = atom->nlocal;
int m = 0;
m = n = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
buf[m++] = tag[i];
@ -151,57 +101,18 @@ int DumpXYZ::pack()
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
if (ids) ids[n++] = tag[i];
}
return m;
}
/* ---------------------------------------------------------------------- */
void DumpXYZ::write_data(int n, double *mybuf)
{
// for group = all, spread buf atom coords into global arrays
// if last chunk of atoms in this snapshot, write global arrays to file
if (igroup == 0) {
int j,tag;
int m = 0;
for (int i = 0; i < n; i++) {
tag = static_cast<int> (mybuf[m]) - 1;
types[tag] = static_cast<int> (mybuf[m+1]);
j = 3*tag;
coords[j++] = mybuf[m+2];
coords[j++] = mybuf[m+3];
coords[j] = mybuf[m+4];
m += size_one;
}
ntotal += n;
if (ntotal == natoms) {
write_frame();
ntotal = 0;
}
// for group != all, write unsorted type,x,y,z to file
} else {
int m = 0;
for (int i = 0; i < n; i++) {
fprintf(fp,format,
static_cast<int> (mybuf[m+1]),mybuf[m+2],mybuf[m+3],mybuf[m+4]);
m += size_one;
}
}
}
/* ---------------------------------------------------------------------- */
void DumpXYZ::write_frame()
{
int m = 0;
for (int i = 0; i < natoms; i++) {
fprintf(fp,format,types[i],coords[m],coords[m+1],coords[m+2]);
m += 3;
for (int i = 0; i < n; i++) {
fprintf(fp,format,
static_cast<int> (mybuf[m+1]),mybuf[m+2],mybuf[m+3],mybuf[m+4]);
m += size_one;
}
}

12
src/dump_xyz.h
View File

@ -27,20 +27,14 @@ namespace LAMMPS_NS {
class DumpXYZ : public Dump {
public:
DumpXYZ(class LAMMPS *, int, char**);
~DumpXYZ();
void init();
double memory_usage();
~DumpXYZ() {}
private:
int natoms,ntotal;
int *types;
float *coords;
void init_style();
void write_header(int);
int count();
int pack();
void pack(int *);
void write_data(int, double *);
void write_frame();
};
}

389
src/irregular.cpp
View File

@ -42,7 +42,11 @@ Irregular::Irregular(LAMMPS *lmp) : Pointers(lmp)
procgrid = comm->procgrid;
grid2proc = comm->grid2proc;
// initialize comm buffers
aplan = NULL;
dplan = NULL;
// initialize buffers for atom comm, not used for datum comm
// these can persist for multiple irregular operations
maxsend = BUFMIN;
buf_send = (double *)
@ -56,6 +60,9 @@ Irregular::Irregular(LAMMPS *lmp) : Pointers(lmp)
Irregular::~Irregular()
{
if (aplan) destroy_atom();
if (dplan) destroy_data();
memory->sfree(buf_send);
memory->sfree(buf_recv);
}
@ -64,7 +71,7 @@ Irregular::~Irregular()
communicate atoms to new owning procs via irregular communication
can be used in place of comm->exchange()
unlike exchange(), allows atoms to have moved arbitrarily long distances
first setup irregular comm pattern, then invoke it
sets up irregular plan, invokes it, destroys it
atoms must be remapped to be inside simulation box before this is called
for triclinic: atoms must be in lamda coords (0-1) before this is called
------------------------------------------------------------------------- */
@ -127,11 +134,10 @@ void Irregular::migrate_atoms()
// create irregular communication plan, perform comm, destroy plan
// returned nrecv = size of buffer needed for incoming atoms
int nrecv;
PlanAtom *plan = create_atom(nsendatom,sizes,proclist,&nrecv);
int nrecv = create_atom(nsendatom,sizes,proclist);
if (nrecv > maxrecv) grow_recv(nrecv);
exchange_atom(plan,buf_send,sizes,buf_recv);
destroy_atom(plan);
exchange_atom(buf_send,sizes,buf_recv);
destroy_atom();
delete [] sizes;
delete [] proclist;
@ -147,22 +153,22 @@ void Irregular::migrate_atoms()
}
/* ----------------------------------------------------------------------
create an irregular communication plan for atoms
create a communication plan for atoms
n = # of atoms to send
sizes = # of doubles for each atom
proclist = proc to send each atom to (none to me)
return nrecvsize = total # of doubles I will recv
proclist = proc to send each atom to (not including self)
return total # of doubles I will recv (not including self)
------------------------------------------------------------------------- */
Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes,
int *proclist, int *nrecvsize)
int Irregular::create_atom(int n, int *sizes, int *proclist)
{
int i;
// allocate plan and work vectors
PlanAtom *plan = (struct PlanAtom *)
memory->smalloc(sizeof(PlanAtom),"irregular:plan");
if (aplan) destroy_atom();
aplan = (struct PlanAtom *)
memory->smalloc(sizeof(PlanAtom),"irregular:aplan");
int *list = new int[nprocs];
int *count = new int[nprocs];
@ -221,7 +227,7 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes,
}
}
// num_send = # of datums I send to each proc
// num_send = # of atoms I send to each proc
for (i = 0; i < nsend; i++) num_send[i] = 0;
for (i = 0; i < n; i++) {
@ -229,11 +235,11 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes,
num_send[isend]++;
}
// count = offsets into n-length index_send for each proc I send to
// index_send = list of which datums to send to each proc
// 1st N1 values are datum indices for 1st proc,
// next N2 values are datum indices for 2nd proc, etc
// offset_send = where each datum starts in send buffer
// count = offsets into index_send for each proc I send to
// index_send = list of which atoms to send to each proc
// 1st N1 values are atom indices for 1st proc,
// next N2 values are atom indices for 2nd proc, etc
// offset_send = where each atom starts in send buffer
count[0] = 0;
for (i = 1; i < nsend; i++) count[i] = count[i-1] + num_send[i-1];
@ -246,7 +252,7 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes,
}
// tell receivers how much data I send
// sendmax = largest # of datums I send in a single message
// sendmax = largest # of doubles I send in a single message
int sendmax = 0;
for (i = 0; i < nsend; i++) {
@ -259,15 +265,15 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes,
// length_recv = total size of message each proc sends me
// nrecvsize = total size of data I recv
*nrecvsize = 0;
int nrecvsize = 0;
for (i = 0; i < nrecv; i++) {
MPI_Recv(&length_recv[i],1,MPI_INT,MPI_ANY_SOURCE,0,world,status);
proc_recv[i] = status->MPI_SOURCE;
*nrecvsize += length_recv[i];
nrecvsize += length_recv[i];
}
// barrier to insure all MPI_ANY_SOURCE messages are received
// else another proc could proceed to irregular_perform() and send to me
// else another proc could proceed to exchange_atom() and send to me
MPI_Barrier(world);
@ -276,41 +282,115 @@ Irregular::PlanAtom *Irregular::create_atom(int n, int *sizes,
delete [] count;
delete [] list;
// initialize plan and return it
// initialize plan
plan->nsend = nsend;
plan->nrecv = nrecv;
plan->sendmax = sendmax;
aplan->nsend = nsend;
aplan->nrecv = nrecv;
aplan->sendmax = sendmax;
plan->proc_send = proc_send;
plan->length_send = length_send;
plan->num_send = num_send;
plan->index_send = index_send;
plan->offset_send = offset_send;
plan->proc_recv = proc_recv;
plan->length_recv = length_recv;
plan->request = request;
plan->status = status;
aplan->proc_send = proc_send;
aplan->length_send = length_send;
aplan->num_send = num_send;
aplan->index_send = index_send;
aplan->offset_send = offset_send;
aplan->proc_recv = proc_recv;
aplan->length_recv = length_recv;
return plan;
aplan->request = request;
aplan->status = status;
return nrecvsize;
}
/* ----------------------------------------------------------------------
create an irregular communication plan for datums
n = # of datums to send
proclist = proc to send each datum to (none to me)
return nrecvsize = total # of datums I will recv
communicate atoms via PlanAtom
sendbuf = list of atoms to send
sizes = # of doubles for each atom
recvbuf = received atoms
------------------------------------------------------------------------- */
Irregular::PlanData *Irregular::create_data(int n, int *proclist,
int *nrecvsize)
void Irregular::exchange_atom(double *sendbuf, int *sizes, double *recvbuf)
{
int i;
int i,m,n,offset,num_send;
// post all receives
offset = 0;
for (int irecv = 0; irecv < aplan->nrecv; irecv++) {
MPI_Irecv(&recvbuf[offset],aplan->length_recv[irecv],MPI_DOUBLE,
aplan->proc_recv[irecv],0,world,&aplan->request[irecv]);
offset += aplan->length_recv[irecv];
}
// allocate buf for largest send
double *buf = (double *) memory->smalloc(aplan->sendmax*sizeof(double),
"irregular:buf");
// send each message
// pack buf with list of atoms
// m = index of atom in sendbuf
int *index_send = aplan->index_send;
int nsend = aplan->nsend;
n = 0;
for (int isend = 0; isend < nsend; isend++) {
offset = 0;
num_send = aplan->num_send[isend];
for (i = 0; i < num_send; i++) {
m = index_send[n++];
memcpy(&buf[offset],&sendbuf[aplan->offset_send[m]],
sizes[m]*sizeof(double));
offset += sizes[m];
}
MPI_Send(buf,aplan->length_send[isend],MPI_DOUBLE,
aplan->proc_send[isend],0,world);
}
// free temporary send buffer
memory->sfree(buf);
// wait on all incoming messages
if (aplan->nrecv) MPI_Waitall(aplan->nrecv,aplan->request,aplan->status);
}
/* ----------------------------------------------------------------------
destroy communication plan for atoms
------------------------------------------------------------------------- */
void Irregular::destroy_atom()
{
delete [] aplan->proc_send;
delete [] aplan->length_send;
delete [] aplan->num_send;
delete [] aplan->index_send;
delete [] aplan->offset_send;
delete [] aplan->proc_recv;
delete [] aplan->length_recv;
delete [] aplan->request;
delete [] aplan->status;
memory->sfree(aplan);
aplan = NULL;
}
/* ----------------------------------------------------------------------
create a communication plan for datums
n = # of datums to send
proclist = proc to send each datum to (including self)
return total # of datums I will recv (including self)
------------------------------------------------------------------------- */
int Irregular::create_data(int n, int *proclist)
{
int i,m;
// allocate plan and work vectors
PlanData *plan = (struct PlanData *)
memory->smalloc(sizeof(PlanData),"irregular:plan");
dplan = (struct PlanData *)
memory->smalloc(sizeof(PlanData),"irregular:dplan");
int *list = new int[nprocs];
int *count = new int[nprocs];
@ -324,6 +404,7 @@ Irregular::PlanData *Irregular::create_data(int n, int *proclist,
int nrecv;
MPI_Reduce_scatter(list,&nrecv,count,MPI_INT,MPI_SUM,world);
if (list[me]) nrecv--;
// allocate receive arrays
@ -340,42 +421,56 @@ Irregular::PlanData *Irregular::create_data(int n, int *proclist,
int nsend = 0;
for (i = 0; i < nprocs; i++)
if (list[i]) nsend++;
if (list[me]) nsend--;
// allocate send arrays
// allocate send and self arrays
int *proc_send = new int[nsend];
int *num_send = new int[nsend];
int *index_send = new int[n];
int *index_send = new int[n-list[me]];
int *index_self = new int[list[me]];
// proc_send = procs I send to
// num_send = # of datums I send to each proc
// num_self = # of datums I copy to self
// to balance pattern of send messages:
// each proc begins with iproc > me, continues until iproc = me
// reset list to store which send message each proc corresponds to
int num_self;
int iproc = me;
int isend = 0;
for (i = 0; i < nprocs; i++) {
iproc++;
if (iproc == nprocs) iproc = 0;
if (list[iproc] > 0) {
if (iproc == me) num_self = list[iproc];
else if (list[iproc] > 0) {
proc_send[isend] = iproc;
num_send[isend] = list[iproc];
list[iproc] = isend;
isend++;
}
}
list[me] = 0;
// count = offsets into n-length index_send for each proc I send to
// count = offsets into index_send for each proc I send to
// m = ptr into index_self
// index_send = list of which datums to send to each proc
// 1st N1 values are datum indices for 1st proc,
// next N2 values are datum indices for 2nd proc, etc
// offset_send = where each datum starts in send buffer
count[0] = 0;
for (i = 1; i < nsend; i++) count[i] = count[i-1] + num_send[i-1];
m = 0;
for (i = 0; i < n; i++) {
isend = list[proclist[i]];
index_send[count[isend]++] = i;
iproc = proclist[i];
if (iproc == me) index_self[m++] = i;
else {
isend = list[iproc];
index_send[count[isend]++] = i;
}
}
// tell receivers how much data I send
@ -392,15 +487,16 @@ Irregular::PlanData *Irregular::create_data(int n, int *proclist,
// num_recv = total size of message each proc sends me
// nrecvsize = total size of data I recv
*nrecvsize = 0;
int nrecvsize = 0;
for (i = 0; i < nrecv; i++) {
MPI_Recv(&num_recv[i],1,MPI_INT,MPI_ANY_SOURCE,0,world,status);
proc_recv[i] = status->MPI_SOURCE;
*nrecvsize += num_recv[i];
nrecvsize += num_recv[i];
}
nrecvsize += num_self;
// barrier to insure all MPI_ANY_SOURCE messages are received
// else another proc could proceed to irregular_perform() and send to me
// else another proc could proceed to exchange_data() and send to me
MPI_Barrier(world);
@ -411,166 +507,105 @@ Irregular::PlanData *Irregular::create_data(int n, int *proclist,
// initialize plan and return it
plan->nsend = nsend;
plan->nrecv = nrecv;
plan->sendmax = sendmax;
dplan->nsend = nsend;
dplan->nrecv = nrecv;
dplan->sendmax = sendmax;
plan->proc_send = proc_send;
plan->num_send = num_send;
plan->index_send = index_send;
plan->proc_recv = proc_recv;
plan->num_recv = num_recv;
plan->request = request;
plan->status = status;
dplan->proc_send = proc_send;
dplan->num_send = num_send;
dplan->index_send = index_send;
dplan->proc_recv = proc_recv;
dplan->num_recv = num_recv;
dplan->num_self = num_self;
dplan->index_self = index_self;
return plan;
dplan->request = request;
dplan->status = status;
return nrecvsize;
}
/* ----------------------------------------------------------------------
perform irregular communication of atoms
plan = previouly computed PlanAtom via create_atom()
sendbuf = list of atoms to send
sizes = # of doubles for each atom
recvbuf = received atoms
------------------------------------------------------------------------- */
void Irregular::exchange_atom(PlanAtom *plan, double *sendbuf, int *sizes,
double *recvbuf)
{
int i,m,n,offset;
// post all receives
offset = 0;
for (int irecv = 0; irecv < plan->nrecv; irecv++) {
MPI_Irecv(&recvbuf[offset],plan->length_recv[irecv],MPI_DOUBLE,
plan->proc_recv[irecv],0,world,&plan->request[irecv]);
offset += plan->length_recv[irecv];
}
// allocate buf for largest send
double *buf = (double *) memory->smalloc(plan->sendmax*sizeof(double),
"irregular:buf");
// send each message
// pack buf with list of datums (datum = one atom)
// m = index of datum in sendbuf
n = 0;
for (int isend = 0; isend < plan->nsend; isend++) {
offset = 0;
for (i = 0; i < plan->num_send[isend]; i++) {
m = plan->index_send[n++];
memcpy(&buf[offset],&sendbuf[plan->offset_send[m]],
sizes[m]*sizeof(double));
offset += sizes[m];
}
MPI_Send(buf,plan->length_send[isend],MPI_DOUBLE,
plan->proc_send[isend],0,world);
}
// free temporary send buffer
memory->sfree(buf);
// wait on all incoming messages
if (plan->nrecv) MPI_Waitall(plan->nrecv,plan->request,plan->status);
}
/* ----------------------------------------------------------------------
perform irregular communication of datums
plan = previouly computed PlanData via create_data()
communicate datums via PlanData
sendbuf = list of datums to send
nbytes = size of each datum
recvbuf = received datums
recvbuf = received datums (including copied from me)
------------------------------------------------------------------------- */
void Irregular::exchange_data(PlanData *plan, char *sendbuf, int nbytes,
char *recvbuf)
void Irregular::exchange_data(char *sendbuf, int nbytes, char *recvbuf)
{
int i,m,n,offset;
int i,m,n,offset,num_send;
// post all receives
// post all receives, starting after self copies
offset = 0;
for (int irecv = 0; irecv < plan->nrecv; irecv++) {
MPI_Irecv(&recvbuf[offset],nbytes*plan->num_recv[irecv],MPI_CHAR,
plan->proc_recv[irecv],0,world,&plan->request[irecv]);
offset += nbytes*plan->num_recv[irecv];
offset = dplan->num_self*nbytes;
for (int irecv = 0; irecv < dplan->nrecv; irecv++) {
MPI_Irecv(&recvbuf[offset],dplan->num_recv[irecv]*nbytes,MPI_CHAR,
dplan->proc_recv[irecv],0,world,&dplan->request[irecv]);
offset += dplan->num_recv[irecv]*nbytes;
}
// allocate buf for largest send
char *buf = (char *) memory->smalloc(plan->sendmax*nbytes,"irregular:buf");
char *buf = (char *) memory->smalloc(dplan->sendmax*nbytes,"irregular:buf");
// send each message
// pack buf with list of datums (datum = one atom)
// pack buf with list of datums
// m = index of datum in sendbuf
int *index_send = dplan->index_send;
int nsend = dplan->nsend;
n = 0;
for (int isend = 0; isend < plan->nsend; isend++) {
for (i = 0; i < plan->num_send[isend]; i++) {
m = plan->index_send[n++];
for (int isend = 0; isend < nsend; isend++) {
num_send = dplan->num_send[isend];
for (i = 0; i < num_send; i++) {
m = index_send[n++];
memcpy(&buf[i*nbytes],&sendbuf[m*nbytes],nbytes);
}
MPI_Send(buf,nbytes*plan->num_send[isend],MPI_CHAR,
plan->proc_send[isend],0,world);
MPI_Send(buf,dplan->num_send[isend]*nbytes,MPI_CHAR,
dplan->proc_send[isend],0,world);
}
// free temporary send buffer
memory->sfree(buf);
// copy datums owned by self
// copy datums to self, put at beginning of recvbuf
for (i = 0; i < plan->num_self; i++) {
m = plan->index_self[i++];
int *index_self = dplan->index_self;
int num_self = dplan->num_self;
for (i = 0; i < num_self; i++) {
m = index_self[i];
memcpy(&recvbuf[i*nbytes],&sendbuf[m*nbytes],nbytes);
}
// wait on all incoming messages
if (plan->nrecv) MPI_Waitall(plan->nrecv,plan->request,plan->status);
if (dplan->nrecv) MPI_Waitall(dplan->nrecv,dplan->request,dplan->status);
}
/* ----------------------------------------------------------------------
destroy an irregular communication plan for atoms
destroy communication plan for datums
------------------------------------------------------------------------- */
void Irregular::destroy_atom(PlanAtom *plan)
void Irregular::destroy_data()
{
delete [] plan->proc_send;
delete [] plan->length_send;
delete [] plan->num_send;
delete [] plan->index_send;
delete [] plan->offset_send;
delete [] plan->proc_recv;
delete [] plan->length_recv;
delete [] plan->request;
delete [] plan->status;
memory->sfree(plan);
delete [] dplan->proc_send;
delete [] dplan->num_send;
delete [] dplan->index_send;
delete [] dplan->proc_recv;
delete [] dplan->num_recv;
delete [] dplan->index_self;
delete [] dplan->request;
delete [] dplan->status;
memory->sfree(dplan);
dplan = NULL;
}
/* ----------------------------------------------------------------------
destroy an irregular communication plan for datums
------------------------------------------------------------------------- */
void Irregular::destroy_data(PlanData *plan)
{
delete [] plan->proc_send;
delete [] plan->num_send;
delete [] plan->index_send;
delete [] plan->proc_recv;
delete [] plan->num_recv;
delete [] plan->request;
delete [] plan->status;
memory->sfree(plan);
}
/* ----------------------------------------------------------------------
determine which proc owns atom with x coord
determine which proc owns atom with coord x[3]
x will be in box (orthogonal) or lamda coords (triclinic)
------------------------------------------------------------------------- */
@ -634,3 +669,13 @@ void Irregular::grow_recv(int n)
"comm:buf_recv");
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
------------------------------------------------------------------------- */
double Irregular::memory_usage()
{
double bytes = maxsend * sizeof(double);
bytes += maxrecv * sizeof(double);
return bytes;
}

61
src/irregular.h
View File

@ -20,21 +20,45 @@ namespace LAMMPS_NS {
class Irregular : protected Pointers {
public:
struct PlanAtom { // plan for irregular communication of atoms
Irregular(class LAMMPS *);
~Irregular();
void migrate_atoms();
int create_data(int, int *);
void exchange_data(char *, int, char *);
void destroy_data();
double memory_usage();
private:
int me,nprocs;
int triclinic;
int map_style;
int *procgrid;
int ***grid2proc;
int maxsend,maxrecv; // size of buffers in # of doubles
double *buf_send,*buf_recv;
// plan for irregular communication of atoms
// no params refer to atoms copied to self
struct PlanAtom {
int nsend; // # of messages to send
int nrecv; // # of messages to recv
int sendmax; // # of doubles in largest send message
int *proc_send; // procs to send to
int *length_send; // # of doubles to send to each proc
int *num_send; // # of datums to send to each proc
int *index_send; // list of which datums to send to each proc
int *offset_send; // where each datum starts in send buffer
int *num_send; // # of atoms to send to each proc
int *index_send; // list of which atoms to send to each proc
int *offset_send; // where each atom starts in send buffer
int *proc_recv; // procs to recv from
int *length_recv; // # of doubles to recv from each proc
MPI_Request *request; // MPI requests for posted recvs
MPI_Status *status; // MPI statuses for WaitAll
};
// plan for irregular communication of datums
// only 2 self params refer to atoms copied to self
struct PlanData { // plan for irregular communication of data
int nsend; // # of messages to send
int nrecv; // # of messages to recv
@ -43,32 +67,19 @@ class Irregular : protected Pointers {
int *num_send; // # of datums to send to each proc
int *index_send; // list of which datums to send to each proc
int *proc_recv; // procs to recv from
int *num_recv; // # of doubles to recv from each proc
int num_self;
int *index_self;
int *num_recv; // # of datums to recv from each proc
int num_self; // # of datums to copy to self
int *index_self; // list of which datums to copy to self
MPI_Request *request; // MPI requests for posted recvs
MPI_Status *status; // MPI statuses for WaitAll
};
Irregular(class LAMMPS *);
~Irregular();
void migrate_atoms();
struct PlanData *create_data(int, int *, int *);
void exchange_data(PlanData *, char *, int, char *);
void destroy_data(PlanData *);
PlanAtom *aplan;
PlanData *dplan;
private:
int me,nprocs;
int triclinic;
int map_style;
int *procgrid;
int ***grid2proc;
int maxsend,maxrecv;
double *buf_send,*buf_recv;
struct PlanAtom *create_atom(int, int *, int *, int *);
void exchange_atom(PlanAtom *, double *, int *, double *);
void destroy_atom(PlanAtom *);
int create_atom(int, int *, int *);
void exchange_atom(double *, int *, double *);
void destroy_atom();
int coord2proc(double *);
void grow_send(int,int); // reallocate send buffer