diff --git a/src/XTC/dump_xtc.h b/src/XTC/dump_xtc.h
index 8d7d89df4..081891b22 100644
--- a/src/XTC/dump_xtc.h
+++ b/src/XTC/dump_xtc.h
@@ -1,107 +1,107 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef DUMP_CLASS
DumpStyle(xtc,DumpXTC)
#else
#ifndef LMP_DUMP_XTC_H
#define LMP_DUMP_XTC_H
#include "dump.h"
#ifdef LAMMPS_XDR
#include "xdr_compat.h"
#else
#include "rpc/rpc.h"
#include "rpc/xdr.h"
#endif
namespace LAMMPS_NS {
class DumpXTC : public Dump {
public:
DumpXTC(class LAMMPS *, int, char**);
- ~DumpXTC();
+ virtual ~DumpXTC();
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(bigint);
void pack(int *);
void write_data(int, double *);
bigint memory_usage();
void write_frame();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Invalid dump xtc filename
Filenames used with the dump xtc style cannot be binary or compressed
or cause multiple files to be written.
E: Too many atoms for dump xtc
The system size must fit in a 32-bit integer to use this dump
style.
E: Dump xtc requires sorting by atom ID
Use the dump_modify sort command to enable this.
E: Cannot set dump_modify flush for dump xtc
Self-explanatory.
E: Cannot use variable every setting for dump xtc
The format of this file requires snapshots at regular intervals.
E: Cannot change dump_modify every for dump xtc
The frequency of writing dump xtc snapshots cannot be changed.
E: Cannot open dump file
The output file for the dump command cannot be opened. Check that the
path and name are correct.
E: Too big a timestep for dump xtc
The timestep must fit in a 32-bit integer to use this dump style.
*/
diff --git a/src/dump.cpp b/src/dump.cpp
index 551cde0ad..9b1b10fe9 100644
--- a/src/dump.cpp
+++ b/src/dump.cpp
@@ -1,810 +1,814 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "lmptype.h"
#include "mpi.h"
#include "stdlib.h"
#include "string.h"
#include "stdio.h"
#include "dump.h"
#include "atom.h"
#include "irregular.h"
#include "update.h"
#include "domain.h"
#include "group.h"
#include "output.h"
#include "memory.h"
#include "error.h"
#include "force.h"
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
enum{ASCEND,DESCEND};
/* ---------------------------------------------------------------------- */
Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
{
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
int n = strlen(arg[0]) + 1;
id = new char[n];
strcpy(id,arg[0]);
igroup = group->find(arg[1]);
groupbit = group->bitmask[igroup];
n = strlen(arg[2]) + 1;
style = new char[n];
strcpy(style,arg[2]);
n = strlen(arg[4]) + 1;
filename = new char[n];
strcpy(filename,arg[4]);
comm_forward = comm_reverse = 0;
first_flag = 0;
flush_flag = 1;
format = NULL;
format_user = NULL;
format_default = NULL;
clearstep = 0;
sort_flag = 0;
append_flag = 0;
padflag = 0;
maxbuf = maxids = maxsort = maxproc = 0;
buf = bufsort = NULL;
ids = idsort = index = proclist = NULL;
irregular = NULL;
// parse filename for special syntax
// if contains '%', write one file per proc and replace % with proc-ID
// if contains '*', write one file per timestep and replace * with timestep
// check file suffixes
// if ends in .bin = binary file
// else if ends in .gz = gzipped text file
// else ASCII text file
fp = NULL;
singlefile_opened = 0;
compressed = 0;
binary = 0;
multifile = 0;
multiproc = 0;
nclusterprocs = nprocs;
filewriter = 0;
if (me == 0) filewriter = 1;
fileproc = 0;
multiname = NULL;
char *ptr;
if (ptr = strchr(filename,'%')) {
multiproc = 1;
nclusterprocs = 1;
filewriter = 1;
fileproc = me;
MPI_Comm_split(world,me,0,&clustercomm);
multiname = new char[strlen(filename) + 16];
*ptr = '\0';
sprintf(multiname,"%s%d%s",filename,me,ptr+1);
*ptr = '%';
}
if (strchr(filename,'*')) multifile = 1;
char *suffix = filename + strlen(filename) - strlen(".bin");
if (suffix > filename && strcmp(suffix,".bin") == 0) binary = 1;
suffix = filename + strlen(filename) - strlen(".gz");
if (suffix > filename && strcmp(suffix,".gz") == 0) compressed = 1;
}
/* ---------------------------------------------------------------------- */
Dump::~Dump()
{
delete [] id;
delete [] style;
delete [] filename;
delete [] multiname;
delete [] format;
delete [] format_default;
delete [] format_user;
memory->destroy(buf);
memory->destroy(bufsort);
memory->destroy(ids);
memory->destroy(idsort);
memory->destroy(index);
memory->destroy(proclist);
delete irregular;
if (multiproc) MPI_Comm_free(&clustercomm);
// XTC style sets fp to NULL since it closes file in its destructor
if (multifile == 0 && fp != NULL) {
if (compressed) {
if (filewriter) pclose(fp);
} else {
if (filewriter) fclose(fp);
}
}
}
/* ---------------------------------------------------------------------- */
void Dump::init()
{
init_style();
if (!sort_flag) {
memory->destroy(bufsort);
memory->destroy(ids);
memory->destroy(idsort);
memory->destroy(index);
memory->destroy(proclist);
delete irregular;
maxids = maxsort = maxproc = 0;
bufsort = NULL;
ids = idsort = index = proclist = NULL;
irregular = NULL;
}
if (sort_flag) {
if (multiproc > 1)
error->all(FLERR,
"Cannot dump sort when multiple procs write the dump file");
if (sortcol == 0 && atom->tag_enable == 0)
error->all(FLERR,"Cannot dump sort on atom IDs with no atom IDs defined");
if (sortcol && sortcol > size_one)
error->all(FLERR,"Dump sort column is invalid");
if (nprocs > 1 && irregular == NULL)
irregular = new Irregular(lmp);
bigint size = group->count(igroup);
if (size > MAXSMALLINT) error->all(FLERR,"Too many atoms to dump sort");
// set reorderflag = 1 if can simply reorder local atoms rather than sort
// criteria: sorting by ID, atom IDs are consecutive from 1 to Natoms
// min/max IDs of group match size of group
// compute ntotal_reorder, nme_reorder, idlo/idhi to test against later
reorderflag = 0;
if (sortcol == 0 && atom->tag_consecutive()) {
int *tag = atom->tag;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int min = IBIG;
int max = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
min = MIN(min,tag[i]);
max = MAX(max,tag[i]);
}
int minall,maxall;
MPI_Allreduce(&min,&minall,1,MPI_INT,MPI_MIN,world);
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
int isize = static_cast<int> (size);
if (maxall-minall+1 == isize) {
reorderflag = 1;
double range = maxall-minall + EPSILON;
idlo = static_cast<int> (range*me/nprocs + minall);
int idhi = static_cast<int> (range*(me+1)/nprocs + minall);
int lom1 = static_cast<int> ((idlo-1-minall)/range * nprocs);
int lo = static_cast<int> ((idlo-minall)/range * nprocs);
int him1 = static_cast<int> ((idhi-1-minall)/range * nprocs);
int hi = static_cast<int> ((idhi-minall)/range * nprocs);
if (me && me == lom1) idlo--;
else if (me && me != lo) idlo++;
if (me+1 == him1) idhi--;
else if (me+1 != hi) idhi++;
nme_reorder = idhi-idlo;
ntotal_reorder = isize;
}
}
}
}
/* ---------------------------------------------------------------------- */
int Dump::count()
{
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;
}
/* ---------------------------------------------------------------------- */
void Dump::write()
{
// if file per timestep, open new file
if (multifile) openfile();
// simulation box bounds
if (domain->triclinic == 0) {
boxxlo = domain->boxlo[0];
boxxhi = domain->boxhi[0];
boxylo = domain->boxlo[1];
boxyhi = domain->boxhi[1];
boxzlo = domain->boxlo[2];
boxzhi = domain->boxhi[2];
} 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];
boxxy = domain->xy;
boxxz = domain->xz;
boxyz = domain->yz;
}
// nme = # of dump lines this proc contributes to dump
nme = count();
// ntotal = total # of dump lines in snapshot
// nmax = max # of dump lines on any proc
bigint bnme = nme;
MPI_Allreduce(&bnme,&ntotal,1,MPI_LMP_BIGINT,MPI_SUM,world);
int nmax;
if (multiproc != nprocs) MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world);
else nmax = nme;
// write timestep header
// for multiproc,
// nheader = # of lines in this file via Allreduce on clustercomm
bigint nheader = ntotal;
if (multiproc)
MPI_Allreduce(&bnme,&nheader,1,MPI_LMP_BIGINT,MPI_SUM,clustercomm);
if (filewriter) write_header(nheader);
// insure filewriter proc can receive everyone's info
// limit nmax*size_one to int since used as arg in MPI_Rsend() below
// pack my data into buf
// if sorting on IDs also request ID list from pack()
// sort buf as needed
if (nmax > maxbuf) {
if ((bigint) nmax * size_one > MAXSMALLINT)
error->all(FLERR,"Too much per-proc info for dump");
maxbuf = nmax;
memory->destroy(buf);
memory->create(buf,maxbuf*size_one,"dump:buf");
}
if (sort_flag && sortcol == 0 && nmax > maxids) {
maxids = nmax;
memory->destroy(ids);
memory->create(ids,maxids,"dump:ids");
}
if (sort_flag && sortcol == 0) pack(ids);
else pack(NULL);
if (sort_flag) sort();
// filewriter = 1 = this proc writes to file
// ping each proc in my cluster, receive its data, write data to file
// else wait for ping from fileproc, send my data to fileproc
int tmp,nlines;
MPI_Status status;
MPI_Request request;
if (filewriter) {
for (int iproc = 0; iproc < nclusterprocs; iproc++) {
if (iproc) {
MPI_Irecv(buf,maxbuf*size_one,MPI_DOUBLE,me+iproc,0,world,&request);
MPI_Send(&tmp,0,MPI_INT,me+iproc,0,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_DOUBLE,&nlines);
nlines /= size_one;
} else nlines = nme;
write_data(nlines,buf);
}
if (flush_flag) fflush(fp);
} else {
MPI_Recv(&tmp,0,MPI_INT,fileproc,0,world,&status);
MPI_Rsend(buf,nme*size_one,MPI_DOUBLE,fileproc,0,world);
}
// if file per timestep, close file if I am filewriter
if (multifile) {
if (compressed) {
if (filewriter) pclose(fp);
} else {
if (filewriter) fclose(fp);
}
}
}
/* ----------------------------------------------------------------------
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 = filename;
if (multiproc) filecurrent = multiname;
if (multifile) {
char *filestar = filecurrent;
filecurrent = new char[strlen(filestar) + 16];
char *ptr = strchr(filestar,'*');
*ptr = '\0';
if (padflag == 0)
sprintf(filecurrent,"%s" BIGINT_FORMAT "%s",
filestar,update->ntimestep,ptr+1);
else {
char bif[8],pad[16];
strcpy(bif,BIGINT_FORMAT);
sprintf(pad,"%%s%%0%d%s%%s",padflag,&bif[1]);
sprintf(filecurrent,pad,filestar,update->ntimestep,ptr+1);
}
*ptr = '*';
}
// each proc with filewriter = 1 opens a file
if (filewriter) {
if (compressed) {
#ifdef LAMMPS_GZIP
char gzip[128];
sprintf(gzip,"gzip -6 > %s",filecurrent);
+#ifdef _WIN32
+ fp = _popen(gzip,"wb");
+#else
fp = popen(gzip,"w");
+#endif
#else
error->one(FLERR,"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(FLERR,"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->destroy(bufsort);
memory->create(bufsort,maxsort*size_one,"dump:bufsort");
memory->destroy(index);
memory->create(index,maxsort,"dump:index");
if (sortcol == 0) {
memory->destroy(idsort);
memory->create(idsort,maxsort,"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->destroy(proclist);
memory->create(proclist,maxproc,"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 + EPSILON;
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->destroy(bufsort);
memory->create(bufsort,maxsort*size_one,"dump:bufsort");
memory->destroy(index);
memory->create(index,maxsort,"dump:index");
if (sortcol == 0) {
memory->destroy(idsort);
memory->create(idsort,maxsort,"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();
}
// if reorder flag is set & total/per-proc counts match pre-computed values,
// then create index directly from idsort
// else quicksort of index using IDs or buf column as comparator
if (reorderflag) {
if (ntotal != ntotal_reorder) reorderflag = 0;
int flag = 0;
if (nme != nme_reorder) flag = 1;
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall) reorderflag = 0;
if (reorderflag)
for (i = 0; i < nme; i++)
index[idsort[i]-idlo] = i;
}
if (!reorderflag) {
dumpptr = this;
for (i = 0; i < nme; i++) index[i] = i;
if (sortcol == 0) qsort(index,nme,sizeof(int),idcompare);
else if (sortorder == ASCEND) qsort(index,nme,sizeof(int),bufcompare);
else qsort(index,nme,sizeof(int),bufcompare_reverse);
}
// reset buf size and maxbuf to largest of any post-sort nme values
// this insures proc 0 can receive everyone's info
int nmax;
MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world);
if (nmax > maxbuf) {
maxbuf = nmax;
memory->destroy(buf);
memory->create(buf,maxbuf*size_one,"dump:buf");
}
// copy data from bufsort to buf using index
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
called via qsort() in sort() method
is a static method so access data via dumpptr
------------------------------------------------------------------------- */
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 values with size_one stride
called via qsort() in sort() method
is a static method so access data via dumpptr
sort in ASCENDing order
------------------------------------------------------------------------- */
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;
}
/* ----------------------------------------------------------------------
compare two buffer values with size_one stride
called via qsort() in sort() method
is a static method so access data via dumpptr
sort in DESCENDing order
------------------------------------------------------------------------- */
int Dump::bufcompare_reverse(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
------------------------------------------------------------------------- */
void Dump::modify_params(int narg, char **arg)
{
if (narg == 0) error->all(FLERR,"Illegal dump_modify command");
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"append") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) append_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) append_flag = 0;
else error->all(FLERR,"Illegal dump_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"every") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
int idump;
for (idump = 0; idump < output->ndump; idump++)
if (strcmp(id,output->dump[idump]->id) == 0) break;
int n;
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
delete [] output->var_dump[idump];
n = strlen(&arg[iarg+1][2]) + 1;
output->var_dump[idump] = new char[n];
strcpy(output->var_dump[idump],&arg[iarg+1][2]);
n = 0;
} else {
n = force->inumeric(FLERR,arg[iarg+1]);
if (n <= 0) error->all(FLERR,"Illegal dump_modify command");
}
output->every_dump[idump] = n;
iarg += 2;
} else if (strcmp(arg[iarg],"first") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) first_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) first_flag = 0;
else error->all(FLERR,"Illegal dump_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"fileper") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
if (!multiproc)
error->all(FLERR,"Cannot use dump_modify fileper "
"without % in dump file name");
int nper = force->inumeric(FLERR,arg[iarg+1]);
if (nper <= 0) error->all(FLERR,"Illegal dump_modify command");
multiproc = nprocs/nper;
if (nprocs % nper) multiproc++;
fileproc = me/nper * nper;
int fileprocnext = MIN(fileproc+nper,nprocs);
nclusterprocs = fileprocnext - fileproc;
if (me == fileproc) filewriter = 1;
else filewriter = 0;
int icluster = fileproc/nper;
MPI_Comm_free(&clustercomm);
MPI_Comm_split(world,icluster,0,&clustercomm);
delete [] multiname;
multiname = new char[strlen(filename) + 16];
char *ptr = strchr(filename,'%');
*ptr = '\0';
sprintf(multiname,"%s%d%s",filename,icluster,ptr+1);
*ptr = '%';
iarg += 2;
} else if (strcmp(arg[iarg],"flush") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[iarg+1],"yes") == 0) flush_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) flush_flag = 0;
else error->all(FLERR,"Illegal dump_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"format") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
delete [] format_user;
format_user = NULL;
if (strcmp(arg[iarg+1],"none")) {
int n = strlen(arg[iarg+1]) + 1;
format_user = new char[n];
strcpy(format_user,arg[iarg+1]);
}
iarg += 2;
} else if (strcmp(arg[iarg],"nfile") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
if (!multiproc)
error->all(FLERR,"Cannot use dump_modify nfile "
"without % in dump file name");
int nfile = force->inumeric(FLERR,arg[iarg+1]);
if (nfile <= 0) error->all(FLERR,"Illegal dump_modify command");
nfile = MIN(nfile,nprocs);
multiproc = nfile;
int icluster = static_cast<int> ((bigint) me * nfile/nprocs);
fileproc = static_cast<int> ((bigint) icluster * nprocs/nfile);
int fcluster = static_cast<int> ((bigint) fileproc * nfile/nprocs);
if (fcluster < icluster) fileproc++;
int fileprocnext =
static_cast<int> ((bigint) (icluster+1) * nprocs/nfile);
fcluster = static_cast<int> ((bigint) fileprocnext * nfile/nprocs);
if (fcluster < icluster+1) fileprocnext++;
nclusterprocs = fileprocnext - fileproc;
if (me == fileproc) filewriter = 1;
else filewriter = 0;
MPI_Comm_free(&clustercomm);
MPI_Comm_split(world,icluster,0,&clustercomm);
delete [] multiname;
multiname = new char[strlen(filename) + 16];
char *ptr = strchr(filename,'%');
*ptr = '\0';
sprintf(multiname,"%s%d%s",filename,icluster,ptr+1);
*ptr = '%';
iarg += 2;
} else if (strcmp(arg[iarg],"pad") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
padflag = force->inumeric(FLERR,arg[iarg+1]);
if (padflag < 0) error->all(FLERR,"Illegal dump_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"sort") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
if (strcmp(arg[iarg+1],"off") == 0) sort_flag = 0;
else if (strcmp(arg[iarg+1],"id") == 0) {
sort_flag = 1;
sortcol = 0;
sortorder = ASCEND;
} else {
sort_flag = 1;
sortcol = force->inumeric(FLERR,arg[iarg+1]);
sortorder = ASCEND;
if (sortcol == 0) error->all(FLERR,"Illegal dump_modify command");
if (sortcol < 0) {
sortorder = DESCEND;
sortcol = -sortcol;
}
sortcolm1 = sortcol - 1;
}
iarg += 2;
} else {
int n = modify_param(narg-iarg,&arg[iarg]);
if (n == 0) error->all(FLERR,"Illegal dump_modify command");
iarg += n;
}
}
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
------------------------------------------------------------------------- */
bigint Dump::memory_usage()
{
bigint bytes = memory->usage(buf,size_one*maxbuf);
if (sort_flag) {
if (sortcol == 0) bytes += memory->usage(ids,maxids);
bytes += memory->usage(bufsort,size_one*maxsort);
if (sortcol == 0) bytes += memory->usage(idsort,maxsort);
bytes += memory->usage(index,maxsort);
bytes += memory->usage(proclist,maxproc);
if (irregular) bytes += irregular->memory_usage();
}
return bytes;
}
diff --git a/src/dump.h b/src/dump.h
index 12113df54..3f581767e 100644
--- a/src/dump.h
+++ b/src/dump.h
@@ -1,174 +1,175 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LMP_DUMP_H
#define LMP_DUMP_H
#include "mpi.h"
#include "stdio.h"
#include "pointers.h"
namespace LAMMPS_NS {
class Dump : protected Pointers {
public:
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
int comm_forward; // size of forward communication (0 if none)
int comm_reverse; // size of reverse communication (0 if none)
// 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();
virtual void write();
virtual int pack_comm(int, int *, double *, int, int *) {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 *) {}
void modify_params(int, char **);
virtual bigint memory_usage();
protected:
int me,nprocs; // proc info
char *filename; // user-specified file
int compressed; // 1 if dump file is written compressed, 0 no
int binary; // 1 if dump file is written binary, 0 no
int multifile; // 0 = one big file, 1 = one file per timestep
int multiproc; // 0 = proc 0 writes for all,
// else # of procs writing files
int nclusterprocs; // # of procs in my cluster that write to one file
int filewriter; // 1 if this proc writes a file, else 0
int fileproc; // ID of proc in my cluster who writes to file
char *multiname; // dump filename with % converted to cluster ID
MPI_Comm clustercomm; // MPI communicator within my cluster of procs
int header_flag; // 0 = item, 2 = xyz
int flush_flag; // 0 if no flush, 1 if flush every dump
int sort_flag; // 1 if sorted output
int append_flag; // 1 if open file in append mode, 0 if not
int padflag; // timestep padding in filename
int singlefile_opened; // 1 = one big file, already opened, else 0
int sortcol; // 0 to sort on ID, 1-N on columns
int sortcolm1; // sortcol - 1
int sortorder; // ASCEND or DESCEND
char boundstr[9]; // encoding of boundary flags
char *format_default; // default format string
char *format_user; // format string set by user
char *format; // format string for the file write
FILE *fp; // file to write dump to
int size_one; // # of quantities for one atom
int nme; // # of atoms in this dump from me
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;
bigint ntotal; // total # of per-atom lines in snapshot
int reorderflag; // 1 if OK to reorder instead of sort
int ntotal_reorder; // # of atoms that must be in snapshot
int nme_reorder; // # of atoms I must own in snapshot
int idlo; // lowest ID I own when reordering
int maxbuf; // size of buf
double *buf; // memory for atom quantities
int maxids; // size of ids
int maxsort; // size of bufsort, idsort, index
int maxproc; // size of proclist
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(bigint) = 0;
virtual int count();
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 *);
static int bufcompare_reverse(const void *, const void *);
};
}
#endif
/* ERROR/WARNING messages:
E: Cannot dump sort when multiple procs write the dump file
UNDOCUMENTED
E: Cannot dump sort on atom IDs with no atom IDs defined
Self-explanatory.
E: Dump sort column is invalid
Self-explanatory.
E: Too many atoms to dump sort
Cannot sort when running with more than 2^31 atoms.
E: Too much per-proc info for dump
Number of local atoms times number of columns must fit in a 32-bit
integer for dump.
E: Cannot open gzipped file
-LAMMPS is attempting to open a gzipped version of the specified file
-but was unsuccessful. Check that the path and name are correct.
+LAMMPS was compiled without support for reading and writing gzipped
+files through a pipeline to the gzip program with -DLAMMPS_GZIP.
E: Cannot open dump file
-The output file for the dump command cannot be opened. Check that the
-path and name are correct.
+The specified file cannot be opened. Check that the path and name are
+correct. If the file is a compressed file, also check that the gzip
+executable can be found and run.
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Cannot use dump_modify fileper without % in dump file name
UNDOCUMENTED
E: Cannot use dump_modify nfile without % in dump file name
UNDOCUMENTED
*/
diff --git a/src/dump_image.cpp b/src/dump_image.cpp
index a0063ff30..51c77158e 100644
--- a/src/dump_image.cpp
+++ b/src/dump_image.cpp
@@ -1,1068 +1,1070 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "ctype.h"
#include "stdlib.h"
#include "string.h"
#include "dump_image.h"
#include "image.h"
#include "atom.h"
#include "domain.h"
#include "group.h"
#include "force.h"
#include "comm.h"
#include "input.h"
#include "variable.h"
#include "math_const.h"
#include "error.h"
#include "memory.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define BIG 1.0e20
-enum{PPM,JPG,PNG};
enum{NUMERIC,ATOM,TYPE,ELEMENT,ATTRIBUTE};
enum{STATIC,DYNAMIC};
enum{NO,YES};
/* ---------------------------------------------------------------------- */
DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
DumpCustom(lmp, narg, arg)
{
if (binary || multiproc) error->all(FLERR,"Invalid dump image filename");
// force binary flag on to avoid corrupted output on Windows
binary = 1;
// set filetype based on filename suffix
int n = strlen(filename);
if (strlen(filename) > 4 && strcmp(&filename[n-4],".jpg") == 0)
filetype = JPG;
else if (strlen(filename) > 4 && strcmp(&filename[n-4],".JPG") == 0)
filetype = JPG;
else if (strlen(filename) > 5 && strcmp(&filename[n-5],".jpeg") == 0)
filetype = JPG;
else if (strlen(filename) > 5 && strcmp(&filename[n-5],".JPEG") == 0)
filetype = JPG;
else if (strlen(filename) > 4 && strcmp(&filename[n-4],".png") == 0)
filetype = PNG;
else if (strlen(filename) > 4 && strcmp(&filename[n-4],".PNG") == 0)
filetype = PNG;
else filetype = PPM;
#ifndef LAMMPS_JPEG
if (filetype == JPG)
error->all(FLERR,"Support for writing images in JPEG format not included");
#endif
#ifndef LAMMPS_PNG
if (filetype == PNG)
error->all(FLERR,"Support for writing images in PNG format not included");
#endif
// atom color,diameter settings
if (nfield != 2) error->all(FLERR,"Illegal dump image command");
acolor = ATTRIBUTE;
if (strcmp(arg[5],"type") == 0) acolor = TYPE;
else if (strcmp(arg[5],"element") == 0) acolor = ELEMENT;
adiam = ATTRIBUTE;
if (strcmp(arg[6],"type") == 0) adiam = TYPE;
else if (strcmp(arg[6],"element") == 0) adiam = ELEMENT;
// create Image class with single colormap for atoms
// change defaults for 2d
image = new Image(lmp,1);
if (domain->dimension == 2) {
image->theta = 0.0;
image->phi = 0.0;
image->up[0] = 0.0; image->up[1] = 1.0; image->up[2] = 0.0;
}
// set defaults for optional args
atomflag = YES;
if (atom->nbondtypes == 0) bondflag = NO;
else {
bondflag = YES;
bcolor = ATOM;
bdiam = NUMERIC;
bdiamvalue = 0.5;
}
thetastr = phistr = NULL;
cflag = STATIC;
cx = cy = cz = 0.5;
cxstr = cystr = czstr = NULL;
upxstr = upystr = upzstr = NULL;
zoomstr = NULL;
perspstr = NULL;
boxflag = YES;
boxdiam = 0.02;
axesflag = NO;
// parse optional args
int iarg = ioptional;
while (iarg < narg) {
if (strcmp(arg[iarg],"adiam") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
adiam = NUMERIC;
adiamvalue = force->numeric(FLERR,arg[iarg+1]);
if (adiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command");
iarg += 2;
} else if (strcmp(arg[iarg],"atom") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
if (strcmp(arg[iarg+1],"yes") == 0) atomflag = YES;
else if (strcmp(arg[iarg+1],"no") == 0) atomflag = NO;
else error->all(FLERR,"Illegal dump image command");
iarg += 2;
} else if (strcmp(arg[iarg],"bond") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
if (atom->nbondtypes == 0)
error->all(FLERR,"Dump image bond not allowed with no bond types");
bondflag = YES;
if (strcmp(arg[iarg+1],"none") == 0) bondflag = NO;
else if (strcmp(arg[iarg+1],"atom") == 0) bcolor = ATOM;
else if (strcmp(arg[iarg+1],"type") == 0) bcolor = TYPE;
else error->all(FLERR,"Illegal dump image command");
if (!islower(arg[iarg+2][0])) {
bdiam = NUMERIC;
bdiamvalue = force->numeric(FLERR,arg[iarg+2]);
if (bdiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command");
} else if (strcmp(arg[iarg+2],"atom") == 0) bdiam = ATOM;
else if (strcmp(arg[iarg+2],"type") == 0) bdiam = TYPE;
else if (strcmp(arg[iarg+2],"none") == 0) bondflag = NO;
else error->all(FLERR,"Illegal dump image command");
iarg += 3;
} else if (strcmp(arg[iarg],"size") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
int width = force->inumeric(FLERR,arg[iarg+1]);
int height = force->inumeric(FLERR,arg[iarg+2]);
if (width <= 0 || height <= 0)
error->all(FLERR,"Illegal dump image command");
image->width = width;
image->height = height;
iarg += 3;
} else if (strcmp(arg[iarg],"view") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
int n = strlen(&arg[iarg+1][2]) + 1;
thetastr = new char[n];
strcpy(thetastr,&arg[iarg+1][2]);
} else {
double theta = force->numeric(FLERR,arg[iarg+1]);
if (theta < 0.0 || theta > 180.0)
error->all(FLERR,"Invalid dump image theta value");
theta *= MY_PI/180.0;
image->theta = theta;
}
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
int n = strlen(&arg[iarg+2][2]) + 1;
phistr = new char[n];
strcpy(phistr,&arg[iarg+2][2]);
} else {
double phi = force->numeric(FLERR,arg[iarg+2]);
phi *= MY_PI/180.0;
image->phi = phi;
}
iarg += 3;
} else if (strcmp(arg[iarg],"center") == 0) {
if (iarg+5 > narg) error->all(FLERR,"Illegal dump image command");
if (strcmp(arg[iarg+1],"s") == 0) cflag = STATIC;
else if (strcmp(arg[iarg+1],"d") == 0) cflag = DYNAMIC;
else error->all(FLERR,"Illegal dump image command");
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
int n = strlen(&arg[iarg+2][2]) + 1;
cxstr = new char[n];
strcpy(cxstr,&arg[iarg+2][2]);
cflag = DYNAMIC;
} else cx = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
int n = strlen(&arg[iarg+3][2]) + 1;
cystr = new char[n];
strcpy(cystr,&arg[iarg+3][2]);
cflag = DYNAMIC;
} else cy = force->numeric(FLERR,arg[iarg+3]);
if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) {
int n = strlen(&arg[iarg+4][2]) + 1;
czstr = new char[n];
strcpy(czstr,&arg[iarg+4][2]);
cflag = DYNAMIC;
} else cz = force->numeric(FLERR,arg[iarg+4]);
iarg += 5;
} else if (strcmp(arg[iarg],"up") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
int n = strlen(&arg[iarg+1][2]) + 1;
upxstr = new char[n];
strcpy(upxstr,&arg[iarg+1][2]);
} else image->up[0] = force->numeric(FLERR,arg[iarg+1]);
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
int n = strlen(&arg[iarg+2][2]) + 1;
upystr = new char[n];
strcpy(upystr,&arg[iarg+2][2]);
} else image->up[1] = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
int n = strlen(&arg[iarg+3][2]) + 1;
upzstr = new char[n];
strcpy(upzstr,&arg[iarg+3][2]);
} else image->up[2] = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"zoom") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
int n = strlen(&arg[iarg+1][2]) + 1;
zoomstr = new char[n];
strcpy(zoomstr,&arg[iarg+1][2]);
} else {
double zoom = force->numeric(FLERR,arg[iarg+1]);
if (zoom <= 0.0) error->all(FLERR,"Illegal dump image command");
image->zoom = zoom;
}
iarg += 2;
} else if (strcmp(arg[iarg],"persp") == 0) {
error->all(FLERR,"Dump image persp option is not yet supported");
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
int n = strlen(&arg[iarg+1][2]) + 1;
perspstr = new char[n];
strcpy(perspstr,&arg[iarg+1][2]);
} else {
double persp = force->numeric(FLERR,arg[iarg+1]);
if (persp < 0.0) error->all(FLERR,"Illegal dump image command");
image->persp = persp;
}
iarg += 2;
} else if (strcmp(arg[iarg],"box") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
if (strcmp(arg[iarg+1],"yes") == 0) boxflag = YES;
else if (strcmp(arg[iarg+1],"no") == 0) boxflag = NO;
else error->all(FLERR,"Illegal dump image command");
boxdiam = force->numeric(FLERR,arg[iarg+2]);
if (boxdiam < 0.0) error->all(FLERR,"Illegal dump image command");
iarg += 3;
} else if (strcmp(arg[iarg],"axes") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
if (strcmp(arg[iarg+1],"yes") == 0) axesflag = YES;
else if (strcmp(arg[iarg+1],"no") == 0) axesflag = NO;
else error->all(FLERR,"Illegal dump image command");
axeslen = force->numeric(FLERR,arg[iarg+2]);
axesdiam = force->numeric(FLERR,arg[iarg+3]);
if (axeslen < 0.0 || axesdiam < 0.0)
error->all(FLERR,"Illegal dump image command");
iarg += 4;
} else if (strcmp(arg[iarg],"shiny") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
double shiny = force->numeric(FLERR,arg[iarg+1]);
if (shiny < 0.0 || shiny > 1.0)
error->all(FLERR,"Illegal dump image command");
image->shiny = shiny;
iarg += 2;
} else if (strcmp(arg[iarg],"ssao") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
if (strcmp(arg[iarg+1],"yes") == 0) image->ssao = YES;
else if (strcmp(arg[iarg+1],"no") == 0) image->ssao = NO;
else error->all(FLERR,"Illegal dump image command");
int seed = force->inumeric(FLERR,arg[iarg+2]);
if (seed <= 0) error->all(FLERR,"Illegal dump image command");
image->seed = seed;
double ssaoint = force->numeric(FLERR,arg[iarg+3]);
if (ssaoint < 0.0 || ssaoint > 1.0)
error->all(FLERR,"Illegal dump image command");
image->ssaoint = ssaoint;
iarg += 4;
} else error->all(FLERR,"Illegal dump image command");
}
// allocate image buffer now that image size is known
image->buffers();
// communication neede for bonds colored by atoms
if (bondflag) {
if (bcolor == ATOM || bdiam == ATOM) comm_forward = 3;
else comm_forward = 1;
}
// additional defaults for dump_modify options
diamtype = new double[ntypes+1];
diamelement = new double[ntypes+1];
colortype = new double*[ntypes+1];
colorelement = new double*[ntypes+1];
for (int i = 1; i <= ntypes; i++) {
diamtype[i] = 1.0;
if (i % 6 == 1) colortype[i] = image->color2rgb("red");
else if (i % 6 == 2) colortype[i] = image->color2rgb("green");
else if (i % 6 == 3) colortype[i] = image->color2rgb("blue");
else if (i % 6 == 4) colortype[i] = image->color2rgb("yellow");
else if (i % 6 == 5) colortype[i] = image->color2rgb("aqua");
else if (i % 6 == 0) colortype[i] = image->color2rgb("cyan");
}
if (bondflag) {
bdiamtype = new double[atom->nbondtypes+1];
bcolortype = new double*[atom->nbondtypes+1];
for (int i = 1; i <= atom->nbondtypes; i++) {
bdiamtype[i] = 0.5;
if (i % 6 == 1) bcolortype[i] = image->color2rgb("red");
else if (i % 6 == 2) bcolortype[i] = image->color2rgb("green");
else if (i % 6 == 3) bcolortype[i] = image->color2rgb("blue");
else if (i % 6 == 4) bcolortype[i] = image->color2rgb("yellow");
else if (i % 6 == 5) bcolortype[i] = image->color2rgb("aqua");
else if (i % 6 == 0) bcolortype[i] = image->color2rgb("cyan");
}
} else {
bdiamtype = NULL;
bcolortype = NULL;
}
// viewflag = DYNAMIC if any view parameter is dynamic
viewflag = STATIC;
if (thetastr || phistr || cflag == DYNAMIC ||
upxstr || upystr || upzstr || zoomstr || perspstr) viewflag = DYNAMIC;
box_bounds();
if (cflag == STATIC) box_center();
if (viewflag == STATIC) view_params();
// local data
maxbufcopy = 0;
chooseghost = NULL;
bufcopy = NULL;
}
/* ---------------------------------------------------------------------- */
DumpImage::~DumpImage()
{
delete image;
delete [] diamtype;
delete [] diamelement;
delete [] colortype;
delete [] colorelement;
delete [] bdiamtype;
delete [] bcolortype;
memory->destroy(chooseghost);
memory->destroy(bufcopy);
}
/* ---------------------------------------------------------------------- */
void DumpImage::init_style()
{
if (multifile == 0)
error->all(FLERR,"Dump image requires one snapshot per file");
if (sort_flag) error->all(FLERR,"Dump image cannot perform sorting");
DumpCustom::init_style();
// check variables
if (thetastr) {
thetavar = input->variable->find(thetastr);
if (thetavar < 0)
error->all(FLERR,"Variable name for dump image theta does not exist");
if (!input->variable->equalstyle(thetavar))
error->all(FLERR,"Variable for dump image theta is invalid style");
}
if (phistr) {
phivar = input->variable->find(phistr);
if (phivar < 0)
error->all(FLERR,"Variable name for dump image phi does not exist");
if (!input->variable->equalstyle(phivar))
error->all(FLERR,"Variable for dump image phi is invalid style");
}
if (cxstr) {
cxvar = input->variable->find(cxstr);
if (cxvar < 0)
error->all(FLERR,"Variable name for dump image center does not exist");
if (!input->variable->equalstyle(cxvar))
error->all(FLERR,"Variable for dump image center is invalid style");
}
if (cystr) {
cyvar = input->variable->find(cystr);
if (cyvar < 0)
error->all(FLERR,"Variable name for dump image center does not exist");
if (!input->variable->equalstyle(cyvar))
error->all(FLERR,"Variable for dump image center is invalid style");
}
if (czstr) {
czvar = input->variable->find(czstr);
if (czvar < 0)
error->all(FLERR,"Variable name for dump image center does not exist");
if (!input->variable->equalstyle(czvar))
error->all(FLERR,"Variable for dump image center is invalid style");
}
if (upxstr) {
upxvar = input->variable->find(upxstr);
if (upxvar < 0)
error->all(FLERR,"Variable name for dump image center does not exist");
if (!input->variable->equalstyle(upxvar))
error->all(FLERR,"Variable for dump image center is invalid style");
}
if (upystr) {
upyvar = input->variable->find(upystr);
if (upyvar < 0)
error->all(FLERR,"Variable name for dump image center does not exist");
if (!input->variable->equalstyle(upyvar))
error->all(FLERR,"Variable for dump image center is invalid style");
}
if (upzstr) {
upzvar = input->variable->find(upzstr);
if (upzvar < 0)
error->all(FLERR,"Variable name for dump image center does not exist");
if (!input->variable->equalstyle(upzvar))
error->all(FLERR,"Variable for dump image center is invalid style");
}
if (zoomstr) {
zoomvar = input->variable->find(zoomstr);
if (zoomvar < 0)
error->all(FLERR,"Variable name for dump image zoom does not exist");
if (!input->variable->equalstyle(zoomvar))
error->all(FLERR,"Variable for dump image zoom is invalid style");
}
if (perspstr) {
perspvar = input->variable->find(perspstr);
if (perspvar < 0)
error->all(FLERR,"Variable name for dump image persp does not exist");
if (!input->variable->equalstyle(perspvar))
error->all(FLERR,"Variable for dump image persp is invalid style");
}
// set up type -> element mapping
if (atomflag && acolor == ELEMENT) {
for (int i = 1; i <= ntypes; i++) {
colorelement[i] = image->element2color(typenames[i]);
if (colorelement[i] == NULL)
error->all(FLERR,"Invalid dump image element name");
}
}
if (atomflag && adiam == ELEMENT) {
for (int i = 1; i <= ntypes; i++) {
diamelement[i] = image->element2diam(typenames[i]);
if (diamelement[i] == 0.0)
error->all(FLERR,"Invalid dump image element name");
}
}
}
/* ---------------------------------------------------------------------- */
void DumpImage::write()
{
// open new file
openfile();
// reset box center and view parameters if dynamic
box_bounds();
if (cflag == DYNAMIC) box_center();
if (viewflag == DYNAMIC) view_params();
// nme = # of atoms this proc will contribute to dump
nme = count();
if (nme > maxbuf) {
maxbuf = nme;
memory->destroy(buf);
memory->create(buf,maxbuf*size_one,"dump:buf");
}
// pack buf with color & diameter
pack(NULL);
// set minmax color range if using dynamic atom color map
if (acolor == ATTRIBUTE && image->map_dynamic(0)) {
double two[2],twoall[2];
double lo = BIG;
double hi = -BIG;
int m = 0;
for (int i = 0; i < nchoose; i++) {
lo = MIN(lo,buf[m]);
hi = MAX(hi,buf[m]);
m += size_one;
}
two[0] = -lo;
two[1] = hi;
MPI_Allreduce(two,twoall,2,MPI_DOUBLE,MPI_MAX,world);
int flag = image->map_minmax(0,-twoall[0],twoall[1]);
if (flag) error->all(FLERR,"Invalid color map min/max values");
}
// create image on each proc, then merge them
image->clear();
create_image();
image->merge();
// write image file
if (me == 0) {
if (filetype == JPG) image->write_JPG(fp);
else if (filetype == PNG) image->write_PNG(fp);
else image->write_PPM(fp);
- fclose(fp);
+ if (multifile) {
+ fclose(fp);
+ fp = NULL;
+ }
}
}
/* ----------------------------------------------------------------------
simulation box bounds
------------------------------------------------------------------------- */
void DumpImage::box_bounds()
{
if (domain->triclinic == 0) {
boxxlo = domain->boxlo[0];
boxxhi = domain->boxhi[0];
boxylo = domain->boxlo[1];
boxyhi = domain->boxhi[1];
boxzlo = domain->boxlo[2];
boxzhi = domain->boxhi[2];
} 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];
boxxy = domain->xy;
boxxz = domain->xz;
boxyz = domain->yz;
}
}
/* ----------------------------------------------------------------------
reset view parameters
called once from constructor if view is STATIC
called every snapshot from write() if view is DYNAMIC
------------------------------------------------------------------------- */
void DumpImage::box_center()
{
if (cxstr) cx = input->variable->compute_equal(cxvar);
if (cystr) cy = input->variable->compute_equal(cyvar);
if (czstr) cz = input->variable->compute_equal(czvar);
image->xctr = boxxlo + cx*(boxxhi-boxxlo);
image->yctr = boxylo + cy*(boxyhi-boxylo);
image->zctr = boxzlo + cz*(boxzhi-boxzlo);
}
/* ----------------------------------------------------------------------
reset view parameters in Image class
called once from constructor if view is STATIC
called every snapshot from write() if view is DYNAMIC
------------------------------------------------------------------------- */
void DumpImage::view_params()
{
// view direction theta and phi
if (thetastr) {
double theta = input->variable->compute_equal(thetavar);
if (theta < 0.0 || theta > 180.0)
error->all(FLERR,"Invalid dump image theta value");
theta *= MY_PI/180.0;
image->theta = theta;
}
if (phistr) {
double phi = input->variable->compute_equal(phivar);
phi *= MY_PI/180.0;
image->phi = phi;
}
// up vector
if (upxstr) image->up[0] = input->variable->compute_equal(upxvar);
if (upystr) image->up[1] = input->variable->compute_equal(upyvar);
if (upzstr) image->up[2] = input->variable->compute_equal(upzvar);
// zoom and perspective
if (zoomstr) image->zoom = input->variable->compute_equal(zoomvar);
if (image->zoom <= 0.0) error->all(FLERR,"Invalid dump image zoom value");
if (perspstr) image->persp = input->variable->compute_equal(perspvar);
if (image->persp < 0.0) error->all(FLERR,"Invalid dump image persp value");
// remainder of view setup is internal to Image class
image->view_params(boxxlo,boxxhi,boxylo,boxyhi,boxzlo,boxzhi);
}
/* ----------------------------------------------------------------------
create image for atoms on this proc
every pixel has depth
------------------------------------------------------------------------- */
void DumpImage::create_image()
{
int i,j,m,itype,atom1,atom2;
double diameter,delx,dely,delz;
double *color,*color1,*color2;
double xmid[3];
// render my atoms
if (atomflag) {
double **x = atom->x;
m = 0;
for (i = 0; i < nchoose; i++) {
j = clist[i];
if (acolor == TYPE) {
itype = static_cast<int> (buf[m]);
color = colortype[itype];
} else if (acolor == ELEMENT) {
itype = static_cast<int> (buf[m]);
color = colorelement[itype];
} else if (acolor == ATTRIBUTE) {
color = image->map_value2color(0,buf[m]);
}
if (adiam == NUMERIC) {
diameter = adiamvalue;
} else if (adiam == TYPE) {
itype = static_cast<int> (buf[m+1]);
diameter = diamtype[itype];
} else if (adiam == ELEMENT) {
itype = static_cast<int> (buf[m+1]);
diameter = diamelement[itype];
} else if (adiam == ATTRIBUTE) {
diameter = buf[m+1];
}
image->draw_sphere(x[j],color,diameter);
m += size_one;
}
}
// render bonds for my atoms
// both atoms in bond must be selected for bond to be rendered
// if newton_bond is off, only render bond once
// render bond in 2 pieces if crosses periodic boundary
// if bond is deleted (type = 0), do not render
// if bond is turned off (type < 0), still render
if (bondflag) {
double **x = atom->x;
int *tag = atom->tag;
int **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
int *num_bond = atom->num_bond;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = atom->nlocal + atom->nghost;
int newton_bond = force->newton_bond;
// communicate choose flag for ghost atoms to know if they are selected
// if bcolor/bdiam = ATOM, setup bufcopy to comm atom color/diam attributes
if (nall > maxbufcopy) {
maxbufcopy = atom->nmax;
memory->destroy(chooseghost);
memory->create(chooseghost,maxbufcopy,"dump:chooseghost");
if (comm_forward == 3) {
memory->destroy(bufcopy);
memory->create(bufcopy,maxbufcopy,2,"dump:bufcopy");
}
}
for (i = 0; i < nlocal; i++) chooseghost[i] = choose[i];
if (comm_forward == 3) {
for (i = 0; i < nlocal; i++) bufcopy[i][0] = bufcopy[i][1] = 0.0;
m = 0;
for (i = 0; i < nchoose; i++) {
j = clist[i];
bufcopy[j][0] = buf[m];
bufcopy[j][1] = buf[m+1];
m += size_one;
}
}
comm->forward_comm_dump(this);
for (i = 0; i < nchoose; i++) {
atom1 = clist[i];
for (m = 0; m < num_bond[atom1]; m++) {
atom2 = atom->map(bond_atom[atom1][m]);
if (atom2 < 0 || !chooseghost[atom2]) continue;
if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
if (bond_type[atom1][m] == 0) continue;
if (bcolor == ATOM) {
if (acolor == TYPE) {
color1 = colortype[type[atom1]];
color2 = colortype[type[atom2]];
} else if (acolor == ELEMENT) {
color1 = colorelement[type[atom1]];
color2 = colorelement[type[atom2]];
} else if (acolor == ATTRIBUTE) {
color1 = image->map_value2color(0,bufcopy[atom1][0]);
color2 = image->map_value2color(0,bufcopy[atom2][0]);
}
} else if (bcolor == TYPE) {
itype = bond_type[atom1][m];
if (itype < 0) itype = -itype;
color = bcolortype[itype];
}
if (bdiam == NUMERIC) {
diameter = bdiamvalue;
} else if (bdiam == ATOM) {
if (adiam == NUMERIC) {
diameter = adiamvalue;
} else if (adiam == TYPE) {
diameter = MIN(diamtype[type[atom1]],diamtype[type[atom1]]);
} else if (adiam == ELEMENT) {
diameter = MIN(diamelement[type[atom1]],diamelement[type[atom1]]);
} else if (adiam == ATTRIBUTE) {
diameter = MIN(bufcopy[atom1][1],bufcopy[atom2][1]);
}
} else if (bdiam == TYPE) {
itype = bond_type[atom1][m];
if (itype < 0) itype = -itype;
diameter = bdiamtype[itype];
}
// draw cylinder in 2 pieces if bcolor = ATOM
// or bond crosses periodic boundary
delx = x[atom2][0] - x[atom1][0];
dely = x[atom2][1] - x[atom1][1];
delz = x[atom2][2] - x[atom1][2];
if (bcolor == ATOM || domain->minimum_image_check(delx,dely,delz)) {
domain->minimum_image(delx,dely,delz);
xmid[0] = x[atom1][0] + 0.5*delx;
xmid[1] = x[atom1][1] + 0.5*dely;
xmid[2] = x[atom1][2] + 0.5*delz;
if (bcolor == ATOM)
image->draw_cylinder(x[atom1],xmid,color1,diameter,3);
else image->draw_cylinder(x[atom1],xmid,color,diameter,3);
xmid[0] = x[atom2][0] - 0.5*delx;
xmid[1] = x[atom2][1] - 0.5*dely;
xmid[2] = x[atom2][2] - 0.5*delz;
if (bcolor == ATOM)
image->draw_cylinder(xmid,x[atom2],color2,diameter,3);
else image->draw_cylinder(xmid,x[atom2],color,diameter,3);
} else image->draw_cylinder(x[atom1],x[atom2],color,diameter,3);
}
}
}
// render outline of simulation box, orthogonal or triclinic
if (boxflag) {
double diameter = MIN(boxxhi-boxxlo,boxyhi-boxylo);
if (domain->dimension == 3) diameter = MIN(diameter,boxzhi-boxzlo);
diameter *= boxdiam;
double (*boxcorners)[3];
double box[8][3];
if (domain->triclinic == 0) {
box[0][0] = boxxlo; box[0][1] = boxylo; box[0][2] = boxzlo;
box[1][0] = boxxhi; box[1][1] = boxylo; box[1][2] = boxzlo;
box[2][0] = boxxlo; box[2][1] = boxyhi; box[2][2] = boxzlo;
box[3][0] = boxxhi; box[3][1] = boxyhi; box[3][2] = boxzlo;
box[4][0] = boxxlo; box[4][1] = boxylo; box[4][2] = boxzhi;
box[5][0] = boxxhi; box[5][1] = boxylo; box[5][2] = boxzhi;
box[6][0] = boxxlo; box[6][1] = boxyhi; box[6][2] = boxzhi;
box[7][0] = boxxhi; box[7][1] = boxyhi; box[7][2] = boxzhi;
boxcorners = box;
} else {
domain->box_corners();
boxcorners = domain->corners;
}
image->draw_box(boxcorners,diameter);
}
// render XYZ axes in red/green/blue
// offset by 10% of box size and scale by axeslen
if (axesflag) {
double diameter = MIN(boxxhi-boxxlo,boxyhi-boxylo);
if (domain->dimension == 3) diameter = MIN(diameter,boxzhi-boxzlo);
diameter *= axesdiam;
double (*boxcorners)[3];
double axes[4][3];
if (domain->triclinic == 0) {
axes[0][0] = boxxlo; axes[0][1] = boxylo; axes[0][2] = boxzlo;
axes[1][0] = boxxhi; axes[1][1] = boxylo; axes[1][2] = boxzlo;
axes[2][0] = boxxlo; axes[2][1] = boxyhi; axes[2][2] = boxzlo;
axes[3][0] = boxxlo; axes[3][1] = boxylo; axes[3][2] = boxzhi;
} else {
domain->box_corners();
boxcorners = domain->corners;
axes[0][0] = boxcorners[0][0];
axes[0][1] = boxcorners[0][1];
axes[0][2] = boxcorners[0][2];
axes[1][0] = boxcorners[1][0];
axes[1][1] = boxcorners[1][1];
axes[1][2] = boxcorners[1][2];
axes[2][0] = boxcorners[2][0];
axes[2][1] = boxcorners[2][1];
axes[2][2] = boxcorners[2][2];
axes[3][0] = boxcorners[4][0];
axes[3][1] = boxcorners[4][1];
axes[3][2] = boxcorners[4][2];
}
double offset = MAX(boxxhi-boxxlo,boxyhi-boxylo);
if (domain->dimension == 3) offset = MAX(offset,boxzhi-boxzlo);
offset *= 0.1;
axes[0][0] -= offset; axes[0][1] -= offset; axes[0][2] -= offset;
axes[1][0] -= offset; axes[1][1] -= offset; axes[1][2] -= offset;
axes[2][0] -= offset; axes[2][1] -= offset; axes[2][2] -= offset;
axes[3][0] -= offset; axes[3][1] -= offset; axes[3][2] -= offset;
axes[1][0] = axes[0][0] + axeslen*(axes[1][0]-axes[0][0]);
axes[1][1] = axes[0][1] + axeslen*(axes[1][1]-axes[0][1]);
axes[1][2] = axes[0][2] + axeslen*(axes[1][2]-axes[0][2]);
axes[2][0] = axes[0][0] + axeslen*(axes[2][0]-axes[0][0]);
axes[2][1] = axes[0][1] + axeslen*(axes[2][1]-axes[0][1]);
axes[2][2] = axes[0][2] + axeslen*(axes[2][2]-axes[0][2]);
axes[3][0] = axes[0][0] + axeslen*(axes[3][0]-axes[0][0]);
axes[3][1] = axes[0][1] + axeslen*(axes[3][1]-axes[0][1]);
axes[3][2] = axes[0][2] + axeslen*(axes[3][2]-axes[0][2]);
image->draw_axes(axes,diameter);
}
}
/* ---------------------------------------------------------------------- */
int DumpImage::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i,j,m;
m = 0;
if (comm_forward == 1) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = chooseghost[j];
}
} else {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = chooseghost[j];
buf[m++] = bufcopy[j][0];
buf[m++] = bufcopy[j][1];
}
}
return comm_forward;
}
/* ---------------------------------------------------------------------- */
void DumpImage::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
if (comm_forward == 1)
for (i = first; i < last; i++) chooseghost[i] = static_cast<int> (buf[m++]);
else {
for (i = first; i < last; i++) {
chooseghost[i] = static_cast<int> (buf[m++]);
bufcopy[i][0] = buf[m++];
bufcopy[i][1] = buf[m++];
}
}
}
/* ---------------------------------------------------------------------- */
int DumpImage::modify_param(int narg, char **arg)
{
int n = DumpCustom::modify_param(narg,arg);
if (n) return n;
if (strcmp(arg[0],"acolor") == 0) {
if (narg < 3) error->all(FLERR,"Illegal dump_modify command");
int nlo,nhi;
force->bounds(arg[1],atom->ntypes,nlo,nhi);
// ptrs = list of ncount colornames separated by '/'
int ncount = 1;
char *nextptr;
char *ptr = arg[2];
while (nextptr = strchr(ptr,'/')) {
ptr = nextptr + 1;
ncount++;
}
char **ptrs = new char*[ncount+1];
ncount = 0;
ptrs[ncount++] = strtok(arg[2],"/");
while (ptrs[ncount++] = strtok(NULL,"/"));
ncount--;
// assign each of ncount colors in round-robin fashion to types
int m = 0;
for (int i = nlo; i <= nhi; i++) {
colortype[i] = image->color2rgb(ptrs[m%ncount]);
if (colortype[i] == NULL)
error->all(FLERR,"Invalid color in dump_modify command");
m++;
}
delete [] ptrs;
return 3;
}
if (strcmp(arg[0],"adiam") == 0) {
if (narg < 3) error->all(FLERR,"Illegal dump_modify command");
int nlo,nhi;
force->bounds(arg[1],atom->ntypes,nlo,nhi);
double diam = force->numeric(FLERR,arg[2]);
if (diam <= 0.0) error->all(FLERR,"Illegal dump_modify command");
for (int i = nlo; i <= nhi; i++) diamtype[i] = diam;
return 3;
}
if (strcmp(arg[0],"amap") == 0) {
if (narg < 6) error->all(FLERR,"Illegal dump_modify command");
if (strlen(arg[3]) != 2) error->all(FLERR,"Illegal dump_modify command");
int factor = 2;
if (arg[3][0] == 's') factor = 1;
int nentry = force->inumeric(FLERR,arg[5]);
if (nentry < 1) error->all(FLERR,"Illegal dump_modify command");
int n = 6 + factor*nentry;
if (narg < n) error->all(FLERR,"Illegal dump_modify command");
int flag = image->map_reset(0,n-1,&arg[1]);
if (flag) error->all(FLERR,"Illegal dump_modify command");
return n;
}
if (strcmp(arg[0],"bcolor") == 0) {
if (narg < 3) error->all(FLERR,"Illegal dump_modify command");
if (atom->nbondtypes == 0)
error->all(FLERR,"Dump modify bcolor not allowed with no bond types");
int nlo,nhi;
force->bounds(arg[1],atom->nbondtypes,nlo,nhi);
// ptrs = list of ncount colornames separated by '/'
int ncount = 1;
char *nextptr;
char *ptr = arg[2];
while (nextptr = strchr(ptr,'/')) {
ptr = nextptr + 1;
ncount++;
}
char **ptrs = new char*[ncount+1];
ncount = 0;
ptrs[ncount++] = strtok(arg[2],"/");
while (ptrs[ncount++] = strtok(NULL,"/"));
ncount--;
// assign each of ncount colors in round-robin fashion to types
int m = 0;
for (int i = nlo; i <= nhi; i++) {
bcolortype[i] = image->color2rgb(ptrs[m%ncount]);
if (bcolortype[i] == NULL)
error->all(FLERR,"Invalid color in dump_modify command");
m++;
}
delete [] ptrs;
return 3;
}
if (strcmp(arg[0],"bdiam") == 0) {
if (narg < 3) error->all(FLERR,"Illegal dump_modify command");
if (atom->nbondtypes == 0)
error->all(FLERR,"Dump modify bdiam not allowed with no bond types");
int nlo,nhi;
force->bounds(arg[1],atom->ntypes,nlo,nhi);
double diam = force->numeric(FLERR,arg[2]);
if (diam <= 0.0) error->all(FLERR,"Illegal dump_modify command");
for (int i = nlo; i <= nhi; i++) bdiamtype[i] = diam;
return 3;
}
if (strcmp(arg[0],"backcolor") == 0) {
if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
double *color = image->color2rgb(arg[1]);
if (color == NULL) error->all(FLERR,"Invalid color in dump_modify command");
image->background[0] = static_cast<int> (color[0]*255.0);
image->background[1] = static_cast<int> (color[1]*255.0);
image->background[2] = static_cast<int> (color[2]*255.0);
return 2;
}
if (strcmp(arg[0],"boxcolor") == 0) {
if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
image->boxcolor = image->color2rgb(arg[1]);
if (image->boxcolor == NULL)
error->all(FLERR,"Invalid color in dump_modify command");
return 2;
}
if (strcmp(arg[0],"color") == 0) {
if (narg < 5) error->all(FLERR,"Illegal dump_modify command");
int flag = image->addcolor(arg[1],force->numeric(FLERR,arg[2]),force->numeric(FLERR,arg[3]),force->numeric(FLERR,arg[4]));
if (flag) error->all(FLERR,"Illegal dump_modify command");
return 5;
}
return 0;
}
diff --git a/src/dump_image.h b/src/dump_image.h
index 53230b6cb..9c42c584c 100644
--- a/src/dump_image.h
+++ b/src/dump_image.h
@@ -1,184 +1,190 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef DUMP_CLASS
DumpStyle(image,DumpImage)
#else
#ifndef LMP_DUMP_IMAGE_H
#define LMP_DUMP_IMAGE_H
#include "dump_custom.h"
namespace LAMMPS_NS {
class DumpImage : public DumpCustom {
public:
DumpImage(class LAMMPS *, int, char**);
- ~DumpImage();
+ virtual ~DumpImage();
int pack_comm(int, int *, double *, int, int *);
void unpack_comm(int, int, double *);
- private:
+ protected:
int filetype;
+ enum{PPM,JPG,PNG};
+
int acolor,adiam; // what determines color/diam of atoms
double adiamvalue; // atom diameter value
int atomflag,bondflag; // 0/1 for draw atoms,bonds
int bcolor,bdiam; // what determines color/diam of bonds
double bdiamvalue; // bond diameter value
char *thetastr,*phistr; // variables for view theta,phi
int thetavar,phivar; // index to theta,phi vars
int cflag; // static/dynamic box center
double cx,cy,cz; // fractional box center
char *cxstr,*cystr,*czstr; // variables for box center
int cxvar,cyvar,czvar; // index to box center vars
char *upxstr,*upystr,*upzstr; // view up vector variables
int upxvar,upyvar,upzvar; // index to up vector vars
char *zoomstr,*perspstr; // view zoom and perspective variables
int zoomvar,perspvar; // index to zoom,persp vars
int boxflag,axesflag; // 0/1 for draw box and axes
double boxdiam,axeslen,axesdiam; // params for drawing box and axes
int viewflag; // overall view is static or dynamic
double *diamtype,*diamelement,*bdiamtype; // per-type diameters
double **colortype,**colorelement,**bcolortype; // per-type colors
class Image *image; // class that renders each image
int *chooseghost; // extended choose array for comm
double **bufcopy; // buffer for communicating bond/atom info
int maxbufcopy;
- void init_style();
+ virtual void init_style();
int modify_param(int, char **);
void write();
void box_center();
void view_params();
void box_bounds();
void create_image();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Invalid dump image filename
The file produced by dump image cannot be binary and must
be for a single processor.
-E: Cannot dump JPG file
+E: Support for writing images in JPEG format not included
LAMMPS was not built with the -DLAMMPS_JPEG switch in the Makefile.
+E: Support for writing images in PNG format not included
+
+LAMMPS was not built with the -DLAMMPS_PNG switch in the Makefile.
+
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Dump image bond not allowed with no bond types
Self-explanatory.
E: Invalid dump image theta value
Theta must be between 0.0 and 180.0 inclusive.
E: Dump image persp option is not yet supported
Self-explanatory.
E: Dump image requires one snapshot per file
Use a "*" in the filename.
E: Dump image cannot perform sorting
Self-explanatory.
E: Variable name for dump image theta does not exist
Self-explanatory.
E: Variable for dump image theta is invalid style
Must be an equal-style variable.
E: Variable name for dump image phi does not exist
Self-explanatory.
E: Variable for dump image phi is invalid style
Must be an equal-style variable.
E: Variable name for dump image center does not exist
Self-explanatory.
E: Variable for dump image center is invalid style
Must be an equal-style variable.
E: Variable name for dump image zoom does not exist
Self-explanatory.
E: Variable for dump image zoom is invalid style
Must be an equal-style variable.
E: Variable name for dump image persp does not exist
Self-explanatory.
E: Variable for dump image persp is invalid style
Must be an equal-style variable.
E: Invalid dump image element name
The specified element name was not in the standard list of elements.
See the dump_modify doc page.
E: Invalid dump image zoom value
Zoom value must be > 0.0.
E: Invalid dump image persp value
Persp value must be >= 0.0.
E: Invalid color in dump_modify command
The specified color name was not in the list of recognized colors.
See the dump_modify doc page.
E: Dump modify bcolor not allowed with no bond types
Self-explanatory.
E: Dump modify bdiam not allowed with no bond types
Self-explanatory.
*/
diff --git a/src/dump_movie.cpp b/src/dump_movie.cpp
new file mode 100644
index 000000000..6c11f64e1
--- /dev/null
+++ b/src/dump_movie.cpp
@@ -0,0 +1,106 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "dump_movie.h"
+#include "comm.h"
+#include "force.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+DumpMovie::DumpMovie(LAMMPS *lmp, int narg, char **arg) :
+ DumpImage(lmp, narg, arg)
+{
+ if (multiproc || compressed || multifile)
+ error->all(FLERR,"Invalid dump movie filename");
+
+ filetype = PPM;
+ bitrate = 2000;
+ framerate = 24;
+ fp = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpMovie::openfile()
+{
+ char moviecmd[1024];
+
+ if ((comm->me == 0) && (fp == NULL)) {
+
+#ifdef LAMMPS_FFMPEG
+ sprintf(moviecmd,"ffmpeg -v error -y -r %.2f -f image2pipe -c:v ppm -i - "
+ "-r 24.0 -b:v %dk %s ", framerate, bitrate, filename);
+#else
+ error->one(FLERR,"Cannot generate movie file");
+#endif
+
+#if defined(_WIN32)
+ fp = _popen(moviecmd,"wb");
+#else
+ fp = popen(moviecmd,"w");
+#endif
+
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Failed to open FFmpeg pipeline to file %s",filename);
+ error->one(FLERR,str);
+ }
+ }
+}
+/* ---------------------------------------------------------------------- */
+
+void DumpMovie::init_style()
+{
+ // initialize image style circumventing multifile check
+
+ multifile = 1;
+ DumpImage::init_style();
+ multifile = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int DumpMovie::modify_param(int narg, char **arg)
+{
+ int n = DumpImage::modify_param(narg,arg);
+ if (n) return n;
+
+ if (strcmp(arg[0],"bitrate") == 0) {
+ if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
+ bitrate = force->inumeric(FLERR,arg[1]);
+ if (bitrate <= 0.0) error->all(FLERR,"Illegal dump_modify command");
+ return 2;
+ }
+
+ if (strcmp(arg[0],"framerate") == 0) {
+ if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
+ framerate = force->numeric(FLERR,arg[1]);
+ if ((framerate <= 0.1) || (framerate > 24.0))
+ error->all(FLERR,"Illegal dump_modify framerate command");
+ return 2;
+ }
+
+ return 0;
+}
+
diff --git a/src/dump_movie.h b/src/dump_movie.h
new file mode 100644
index 000000000..2cb4d878d
--- /dev/null
+++ b/src/dump_movie.h
@@ -0,0 +1,73 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef DUMP_CLASS
+
+DumpStyle(movie,DumpMovie)
+
+#else
+
+#ifndef LMP_DUMP_MOVIE_H
+#define LMP_DUMP_MOVIE_H
+
+#include "dump_image.h"
+
+namespace LAMMPS_NS {
+
+class DumpMovie : public DumpImage {
+ public:
+ DumpMovie(LAMMPS *, int, char**);
+
+ virtual void openfile();
+ virtual void init_style();
+ virtual int modify_param(int, char **);
+
+ protected:
+ double framerate; // frame rate of animation
+ int bitrate; // bitrate of video file in kbps
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Invalid dump movie filename
+
+The file produced by dump movie cannot be binary or compressed
+and must be a single file for a single processor.
+
+E: Cannot generate movie file
+
+LAMMPS was built without the -DLAMMPS_FFMPEG switch in the Makefile
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: pipe:: Input/output error
+
+Harmless. This happens when the pipeline to FFmpeg is closed and no
+more image data is sent to be appended to the movie. FFmpeg will
+simply terminate and close the movie file.
+
+E: Failed to open FFmpeg pipeline to file %s
+
+The specified file cannot be opened. Check that the path and name are
+correct and writable and that the FFmpeg executable can be found and run.
+
+*/
diff --git a/src/fix_tmd.cpp b/src/fix_tmd.cpp
index e8ab940d6..11c8c7920 100644
--- a/src/fix_tmd.cpp
+++ b/src/fix_tmd.cpp
@@ -1,541 +1,547 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Paul Crozier (SNL)
Christian Burisch (Bochum Univeristy, Germany)
------------------------------------------------------------------------- */
#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "fix_tmd.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "group.h"
#include "respa.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
#define CHUNK 1000
#define MAXLINE 256
/* ---------------------------------------------------------------------- */
FixTMD::FixTMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (narg < 6) error->all(FLERR,"Illegal fix tmd command");
rho_stop = force->numeric(FLERR,arg[3]);
nfileevery = force->inumeric(FLERR,arg[5]);
if (rho_stop < 0 || nfileevery < 0)
error->all(FLERR,"Illegal fix tmd command");
if (nfileevery && narg != 7) error->all(FLERR,"Illegal fix tmd command");
MPI_Comm_rank(world,&me);
// perform initial allocation of atom-based arrays
// register with Atom class
xf = NULL;
xold = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
// make sure an atom map exists before reading in target coordinates
if (atom->map_style == 0)
error->all(FLERR,"Cannot use fix TMD unless atom map exists");
// read from arg[4] and store coordinates of final target in xf
readfile(arg[4]);
// open arg[6] statistics file and write header
if (nfileevery) {
if (narg != 7) error->all(FLERR,"Illegal fix tmd command");
if (me == 0) {
fp = fopen(arg[6],"w");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open fix tmd file %s",arg[6]);
error->one(FLERR,str);
}
fprintf(fp,"%s %s\n","# Step rho_target rho_old gamma_back",
"gamma_forward lambda work_lambda work_analytical");
}
}
masstotal = group->mass(igroup);
// rho_start = initial rho
// xold = initial x or 0.0 if not in group
int *mask = atom->mask;
int *type = atom->type;
tagint *image = atom->image;
double **x = atom->x;
double *mass = atom->mass;
int nlocal = atom->nlocal;
double dx,dy,dz;
rho_start = 0.0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
domain->unmap(x[i],image[i],xold[i]);
dx = xold[i][0] - xf[i][0];
dy = xold[i][1] - xf[i][1];
dz = xold[i][2] - xf[i][2];
rho_start += mass[type[i]]*(dx*dx + dy*dy + dz*dz);
} else xold[i][0] = xold[i][1] = xold[i][2] = 0.0;
}
double rho_start_total;
MPI_Allreduce(&rho_start,&rho_start_total,1,MPI_DOUBLE,MPI_SUM,world);
rho_start = sqrt(rho_start_total/masstotal);
rho_old = rho_start;
work_lambda = 0.0;
work_analytical = 0.0;
previous_stat = 0;
}
/* ---------------------------------------------------------------------- */
FixTMD::~FixTMD()
{
if (nfileevery && me == 0) fclose(fp);
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,0);
// delete locally stored arrays
memory->destroy(xf);
memory->destroy(xold);
}
/* ---------------------------------------------------------------------- */
int FixTMD::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= INITIAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixTMD::init()
{
// check that no integrator fix comes after a TMD fix
int flag = 0;
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"tmd") == 0) flag = 1;
if (flag && strcmp(modify->fix[i]->style,"nve") == 0) flag = 2;
if (flag && strcmp(modify->fix[i]->style,"nvt") == 0) flag = 2;
if (flag && strcmp(modify->fix[i]->style,"npt") == 0) flag = 2;
if (flag && strcmp(modify->fix[i]->style,"nph") == 0) flag = 2;
}
if (flag == 2) error->all(FLERR,"Fix tmd must come after integration fixes");
// timesteps
dtv = update->dt;
dtf = update->dt * force->ftm2v;
if (strstr(update->integrate_style,"respa"))
step_respa = ((Respa *) update->integrate)->step;
}
/* ---------------------------------------------------------------------- */
void FixTMD::initial_integrate(int vflag)
{
double a,b,c,d,e;
double dx,dy,dz,dxkt,dykt,dzkt;
double dxold,dyold,dzold,xback,yback,zback;
double gamma_forward,gamma_back,gamma_max,lambda;
double kt,fr,kttotal,frtotal,dtfm;
double unwrap[3];
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *mass = atom->mass;
tagint *image = atom->image;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
double rho_target = rho_start + delta * (rho_stop - rho_start);
// compute the Lagrange multiplier
a = b = e = 0.0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dxold = xold[i][0] - xf[i][0];
dyold = xold[i][1] - xf[i][1];
dzold = xold[i][2] - xf[i][2];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xf[i][0];
dy = unwrap[1] - xf[i][1];
dz = unwrap[2] - xf[i][2];
a += mass[type[i]]*(dxold*dxold + dyold*dyold + dzold*dzold);
b += mass[type[i]]*(dx *dxold + dy *dyold + dz *dzold);
e += mass[type[i]]*(dx *dx + dy *dy + dz *dz);
}
}
double abe[3],abetotal[3];
abe[0] = a; abe[1] = b; abe[2] = e;
MPI_Allreduce(abe,abetotal,3,MPI_DOUBLE,MPI_SUM,world);
a = abetotal[0]/masstotal;
b = 2.0*abetotal[1]/masstotal;
e = abetotal[2]/masstotal;
c = e - rho_old*rho_old;
d = b*b - 4*a*c;
if (d < 0) d = 0;
if (b >= 0) gamma_max = (-b - sqrt(d))/(2*a);
else gamma_max = (-b + sqrt(d))/(2*a);
gamma_back = c/(a*gamma_max);
if (a == 0.0) gamma_back = 0;
c = e - rho_target*rho_target;
d = b*b - 4*a*c;
if (d < 0) d = 0;
if (b >= 0) gamma_max = (-b - sqrt(d))/(2*a);
else gamma_max = (-b + sqrt(d))/(2*a);
gamma_forward = c/(a*gamma_max);
if (a == 0.0) gamma_forward = 0;
fr = kt = 0.0;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dxold = xold[i][0] - xf[i][0];
dyold = xold[i][1] - xf[i][1];
dzold = xold[i][2] - xf[i][2];
domain->unmap(x[i],image[i],unwrap);
xback = unwrap[0] + gamma_back*dxold;
yback = unwrap[1] + gamma_back*dyold;
zback = unwrap[2] + gamma_back*dzold;
dxkt = xback - xold[i][0];
dykt = yback - xold[i][1];
dzkt = zback - xold[i][2];
kt += mass[type[i]]*(dxkt*dxkt + dykt*dykt + dzkt*dzkt);
fr += f[i][0]*dxold + f[i][1]*dyold + f[i][2]*dzold;
}
}
double r[2],rtotal[2];
r[0] = fr; r[1] = kt;
MPI_Allreduce(r,rtotal,2,MPI_DOUBLE,MPI_SUM,world);
frtotal = rtotal[0];
kttotal = rtotal[1];
// stat write of mean constraint force based on previous time step constraint
if (nfileevery && me == 0) {
work_analytical +=
(-frtotal - kttotal/dtv/dtf)*(rho_target - rho_old)/rho_old;
lambda = gamma_back*rho_old*masstotal/dtv/dtf;
work_lambda += lambda*(rho_target - rho_old);
if (!(update->ntimestep % nfileevery) &&
(previous_stat != update->ntimestep)) {
fprintf(fp,
BIGINT_FORMAT " %g %g %g %g %g %g %g\n",
update->ntimestep,rho_target,rho_old,
gamma_back,gamma_forward,lambda,work_lambda,work_analytical);
fflush(fp);
previous_stat = update->ntimestep;
}
}
rho_old = rho_target;
// apply the constraint and save constrained positions for next step
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
dxold = xold[i][0] - xf[i][0];
x[i][0] += gamma_forward*dxold;
v[i][0] += gamma_forward*dxold/dtv;
f[i][0] += gamma_forward*dxold/dtv/dtfm;
dyold = xold[i][1] - xf[i][1];
x[i][1] += gamma_forward*dyold;
v[i][1] += gamma_forward*dyold/dtv;
f[i][1] += gamma_forward*dyold/dtv/dtfm;
dzold = xold[i][2] - xf[i][2];
x[i][2] += gamma_forward*dzold;
v[i][2] += gamma_forward*dzold/dtv;
f[i][2] += gamma_forward*dzold/dtv/dtfm;
domain->unmap(x[i],image[i],xold[i]);
}
}
}
/* ---------------------------------------------------------------------- */
void FixTMD::initial_integrate_respa(int vflag, int ilevel, int flag)
{
if (flag) return; // only used by NPT,NPH
dtv = step_respa[ilevel];
dtf = step_respa[ilevel] * force->ftm2v;
if (ilevel == 0) initial_integrate(vflag);
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixTMD::memory_usage()
{
double bytes = 2*atom->nmax*3 * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate atom-based arrays
------------------------------------------------------------------------- */
void FixTMD::grow_arrays(int nmax)
{
memory->grow(xf,nmax,3,"fix_tmd:xf");
memory->grow(xold,nmax,3,"fix_tmd:xold");
}
/* ----------------------------------------------------------------------
copy values within local atom-based arrays
------------------------------------------------------------------------- */
void FixTMD::copy_arrays(int i, int j, int delflag)
{
xf[j][0] = xf[i][0];
xf[j][1] = xf[i][1];
xf[j][2] = xf[i][2];
xold[j][0] = xold[i][0];
xold[j][1] = xold[i][1];
xold[j][2] = xold[i][2];
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for exchange with another proc
------------------------------------------------------------------------- */
int FixTMD::pack_exchange(int i, double *buf)
{
buf[0] = xf[i][0];
buf[1] = xf[i][1];
buf[2] = xf[i][2];
buf[3] = xold[i][0];
buf[4] = xold[i][1];
buf[5] = xold[i][2];
return 6;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based arrays from exchange with another proc
------------------------------------------------------------------------- */
int FixTMD::unpack_exchange(int nlocal, double *buf)
{
xf[nlocal][0] = buf[0];
xf[nlocal][1] = buf[1];
xf[nlocal][2] = buf[2];
xold[nlocal][0] = buf[3];
xold[nlocal][1] = buf[4];
xold[nlocal][2] = buf[5];
return 6;
}
/* ----------------------------------------------------------------------
read target coordinates from file, store with appropriate atom
------------------------------------------------------------------------- */
void FixTMD::readfile(char *file)
{
if (me == 0) {
if (screen) fprintf(screen,"Reading TMD target file %s ...\n",file);
open(file);
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
char *buffer = new char[CHUNK*MAXLINE];
- char *ptr,*next,*bufptr;
+ char *next,*bufptr;
int i,m,nlines,tag,imageflag,ix,iy,iz;
double x,y,z,xprd,yprd,zprd;
int firstline = 1;
int ncount = 0;
char *eof = NULL;
xprd = yprd = zprd = -1.0;
while (!eof) {
if (me == 0) {
m = 0;
for (nlines = 0; nlines < CHUNK; nlines++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) break;
m += strlen(&buffer[m]);
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
MPI_Bcast(&nlines,1,MPI_INT,0,world);
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
bufptr = buffer;
for (i = 0; i < nlines; i++) {
next = strchr(bufptr,'\n');
*next = '\0';
if (firstline) {
if (strstr(bufptr,"xlo xhi")) {
double lo,hi;
sscanf(bufptr,"%lg %lg",&lo,&hi);
xprd = hi - lo;
bufptr = next + 1;
continue;
} else if (strstr(bufptr,"ylo yhi")) {
double lo,hi;
sscanf(bufptr,"%lg %lg",&lo,&hi);
yprd = hi - lo;
bufptr = next + 1;
continue;
} else if (strstr(bufptr,"zlo zhi")) {
double lo,hi;
sscanf(bufptr,"%lg %lg",&lo,&hi);
zprd = hi - lo;
bufptr = next + 1;
continue;
} else if (atom->count_words(bufptr) == 4) {
if (xprd >= 0.0 || yprd >= 0.0 || zprd >= 0.0)
error->all(FLERR,"Incorrect format in TMD target file");
imageflag = 0;
firstline = 0;
} else if (atom->count_words(bufptr) == 7) {
if (xprd < 0.0 || yprd < 0.0 || zprd < 0.0)
error->all(FLERR,"Incorrect format in TMD target file");
imageflag = 1;
firstline = 0;
} else error->all(FLERR,"Incorrect format in TMD target file");
}
if (imageflag)
sscanf(bufptr,"%d %lg %lg %lg %d %d %d",&tag,&x,&y,&z,&ix,&iy,&iz);
else
sscanf(bufptr,"%d %lg %lg %lg",&tag,&x,&y,&z);
m = atom->map(tag);
if (m >= 0 && m < nlocal && mask[m] & groupbit) {
if (imageflag) {
xf[m][0] = x + ix*xprd;
xf[m][1] = y + iy*yprd;
xf[m][2] = z + iz*zprd;
} else {
xf[m][0] = x;
xf[m][1] = y;
xf[m][2] = z;
}
ncount++;
}
bufptr = next + 1;
}
}
// clean up
delete [] buffer;
if (me == 0) {
if (compressed) pclose(fp);
else fclose(fp);
}
// check that all atoms in group were listed in target file
// set xf = 0.0 for atoms not in group
int gcount = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) gcount++;
else xf[i][0] = xf[i][1] = xf[i][2] = 0.0;
int flag = 0;
if (gcount != ncount) flag = 1;
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall) error->all(FLERR,"TMD target file did not list all group atoms");
}
/* ----------------------------------------------------------------------
proc 0 opens TMD data file
test if gzipped
------------------------------------------------------------------------- */
void FixTMD::open(char *file)
{
compressed = 0;
char *suffix = file + strlen(file) - 3;
if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1;
if (!compressed) fp = fopen(file,"r");
else {
#ifdef LAMMPS_GZIP
char gunzip[128];
- sprintf(gunzip,"gunzip -c %s",file);
+ sprintf(gunzip,"gzip -c -d %s",file);
+
+#ifdef _WIN32
+ fp = _popen(gunzip,"rb");
+#else
fp = popen(gunzip,"r");
+#endif
+
#else
error->one(FLERR,"Cannot open gzipped file");
#endif
}
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",file);
error->one(FLERR,str);
}
}
/* ---------------------------------------------------------------------- */
void FixTMD::reset_dt()
{
dtv = update->dt;
dtf = update->dt * force->ftm2v;
}
diff --git a/src/fix_tmd.h b/src/fix_tmd.h
index a6b96fbb1..d00bb303e 100644
--- a/src/fix_tmd.h
+++ b/src/fix_tmd.h
@@ -1,109 +1,110 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(tmd,FixTMD)
#else
#ifndef LMP_FIX_TMD_H
#define LMP_FIX_TMD_H
#include "stdio.h"
#include "fix.h"
namespace LAMMPS_NS {
class FixTMD : public Fix {
public:
FixTMD(class LAMMPS *, int, char **);
~FixTMD();
int setmask();
void init();
void initial_integrate(int);
void initial_integrate_respa(int, int, int);
double memory_usage();
void grow_arrays(int);
void copy_arrays(int, int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
void reset_dt();
private:
int me;
int nfileevery,compressed;
bigint previous_stat;
FILE *fp;
double rho_start,rho_stop,rho_old,masstotal;
double dtv,dtf;
double *step_respa;
double work_lambda,work_analytical;
double **xf,**xold;
void readfile(char *);
void open(char *);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Cannot use fix TMD unless atom map exists
Using this fix requires the ability to lookup an atom index, which is
provided by an atom map. An atom map does not exist (by default) for
non-molecular problems. Using the atom_modify map command will force
an atom map to be created.
E: Cannot open fix tmd file %s
The output file for the fix tmd command cannot be opened. Check that
the path and name are correct.
E: Fix tmd must come after integration fixes
Any fix tmd command must appear in the input script after all time
integration fixes (nve, nvt, npt). See the fix tmd documentation for
details.
E: Incorrect format in TMD target file
Format of file read by fix tmd command is incorrect.
E: TMD target file did not list all group atoms
The target file for the fix tmd command did not list all atoms in the
fix group.
E: Cannot open gzipped file
-LAMMPS is attempting to open a gzipped version of the specified file
-but was unsuccessful. Check that the path and name are correct.
+LAMMPS was compiled without support for reading and writing gzipped
+files through a pipeline to the gzip program with -DLAMMPS_GZIP.
E: Cannot open file %s
The specified file cannot be opened. Check that the path and name are
-correct.
+correct. If the file is a compressed file, also check that the gzip
+executable can be found and run.
*/
diff --git a/src/reader.cpp b/src/reader.cpp
index 79dbb1c04..588b21062 100644
--- a/src/reader.cpp
+++ b/src/reader.cpp
@@ -1,69 +1,75 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "stdio.h"
#include "string.h"
#include "reader.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
Reader::Reader(LAMMPS *lmp) : Pointers(lmp)
{
fp = NULL;
}
/* ----------------------------------------------------------------------
try to open given file
generic version for ASCII files that may be compressed
------------------------------------------------------------------------- */
void Reader::open_file(const char *file)
{
if (fp != NULL) close_file();
compressed = 0;
const char *suffix = file + strlen(file) - 3;
if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1;
if (!compressed) fp = fopen(file,"r");
else {
#ifdef LAMMPS_GZIP
char gunzip[1024];
- sprintf(gunzip,"gunzip -c %s",file);
+ sprintf(gunzip,"gzip -c -d %s",file);
+
+#ifdef _WIN32
+ fp = _popen(gunzip,"rb");
+#else
fp = popen(gunzip,"r");
+#endif
+
#else
error->one(FLERR,"Cannot open gzipped file");
#endif
}
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",file);
error->one(FLERR,str);
}
}
/* ----------------------------------------------------------------------
close current file if open
generic version for ASCII files that may be compressed
------------------------------------------------------------------------- */
void Reader::close_file()
{
if (fp == NULL) return;
if (compressed) pclose(fp);
else fclose(fp);
fp = NULL;
}
diff --git a/src/reader.h b/src/reader.h
index 571016ba7..8f36bf622 100644
--- a/src/reader.h
+++ b/src/reader.h
@@ -1,60 +1,61 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
Contributed by Timothy Sirk
------------------------------------------------------------------------- */
#ifndef LMP_READER_H
#define LMP_READER_H
#include "pointers.h"
namespace LAMMPS_NS {
class Reader : protected Pointers {
public:
Reader(class LAMMPS *);
virtual ~Reader() {}
virtual void settings(int, char**) {};
virtual int read_time(bigint &) = 0;
virtual void skip() = 0;
virtual bigint read_header(double [3][3], int &, int, int, int *, char **,
int, int, int &, int &, int &, int &) = 0;
virtual void read_atoms(int, int, double **) = 0;
virtual void open_file(const char *);
virtual void close_file();
protected:
FILE *fp; // pointer to opened file or pipe
int compressed; // flag for dump file compression
};
}
#endif
/* ERROR/WARNING messages:
E: Cannot open gzipped file
-LAMMPS is attempting to open a gzipped version of the specified file
-but was unsuccessful. Check that the path and name are correct.
+LAMMPS was compiled without support for reading and writing gzipped
+files through a pipeline to the gzip program with -DLAMMPS_GZIP.
E: Cannot open file %s
The specified file cannot be opened. Check that the path and name are
-correct.
+correct. If the file is a compressed file, also check that the gzip
+executable can be found and run.
*/