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

src/STUBS/mpi.cpp

@@ -228,6 +228,25 @@ void MPI_Allreduce(void *sendbuf, void *recvbuf, int count,
 
 /* ---------------------------------------------------------------------- */
 
+/* copy values from data1 to data2 */
+
+void MPI_Reduce(void *sendbuf, void *recvbuf, int count,
+		   MPI_Datatype datatype, MPI_Op op,
+		   int root, MPI_Comm comm)
+{
+  int n;
+  if (datatype == MPI_INT) n = count*sizeof(int);
+  else if (datatype == MPI_FLOAT) n = count*sizeof(float);
+  else if (datatype == MPI_DOUBLE) n = count*sizeof(double);
+  else if (datatype == MPI_CHAR) n = count*sizeof(char);
+  else if (datatype == MPI_BYTE) n = count*sizeof(char);
+  else if (datatype == MPI_DOUBLE_INT) n = count*sizeof(double_int);
+
+  memcpy(recvbuf,sendbuf,n);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void MPI_Scan(void *sendbuf, void *recvbuf, int count,
 	      MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
 {
@@ -316,3 +335,22 @@ void MPI_Gather(void *sendbuf, int sendcount, MPI_Datatype sendtype,
 
   memcpy(recvbuf,sendbuf,n);
 }
+
+/* ---------------------------------------------------------------------- */
+
+/* copy values from data1 to data2 */
+
+void MPI_Gatherv(void *sendbuf, int sendcount, MPI_Datatype sendtype,
+		    void *recvbuf, int *recvcounts, int *displs,
+		    MPI_Datatype recvtype, int root, MPI_Comm comm)
+{
+  int n;
+  if (sendtype == MPI_INT) n = sendcount*sizeof(int);
+  else if (sendtype == MPI_FLOAT) n = sendcount*sizeof(float);
+  else if (sendtype == MPI_DOUBLE) n = sendcount*sizeof(double);
+  else if (sendtype == MPI_CHAR) n = sendcount*sizeof(char);
+  else if (sendtype == MPI_BYTE) n = sendcount*sizeof(char);
+  else if (sendtype == MPI_DOUBLE_INT) n = sendcount*sizeof(double_int);
+
+  memcpy(recvbuf,sendbuf,n);
+}

src/STUBS/mpi.h

@@ -88,6 +88,8 @@ void MPI_Bcast(void *buf, int count, MPI_Datatype datatype,
 	       int root, MPI_Comm comm);
 void MPI_Allreduce(void *sendbuf, void *recvbuf, int count,
 		   MPI_Datatype datatype, MPI_Op op, MPI_Comm comm);
+void MPI_Reduce(void *sendbuf, void *recvbuf, int count,
+		   MPI_Datatype datatype, MPI_Op op, int root, MPI_Comm comm);
 void MPI_Scan(void *sendbuf, void *recvbuf, int count,
 	      MPI_Datatype datatype, MPI_Op op, MPI_Comm comm);
 void MPI_Allgather(void *sendbuf, int sendcount, MPI_Datatype sendtype,
@@ -101,5 +103,8 @@ void MPI_Reduce_scatter(void *sendbuf, void *recvbuf, int *recvcounts,
 void MPI_Gather(void *sendbuf, int sendcount, MPI_Datatype sendtype,
 		void *recvbuf, int recvcount, MPI_Datatype recvtype,
 		int root, MPI_Comm comm);
+void MPI_Gatherv(void *sendbuf, int sendcount, MPI_Datatype sendtype,
+		    void *recvbuf, int *recvcounts, int *displs,
+		    MPI_Datatype recvtype, int root, MPI_Comm comm);
 
 #endif

src/dump.cpp

@@ -52,6 +52,8 @@ Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   format = NULL;
   format_user = NULL;
 
+  sort_flag = 0;
+
   maxbuf = 0;
   buf = NULL;
 
@@ -256,6 +258,12 @@ void Dump::modify_params(int narg, char **arg)
 	if (strcmp(id,output->dump[idump]->id) == 0) break;
       output->dump_every[idump] = n;
       iarg += 2;
+    } else if (strcmp(arg[iarg],"sort") == 0) {
+      if (iarg+2 > narg) error->all("Illegal dump_modify command");
+      if (strcmp(arg[iarg+1],"yes") == 0) sort_flag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) sort_flag = 0;
+      else error->all("Illegal dump_modify command");
+      iarg += 2;
     } else {
       int n = modify_param(narg-iarg,&arg[iarg]);
       if (n == 0) error->all("Illegal dump_modify command");

src/dump.h

@@ -34,6 +34,7 @@ class Dump : protected Pointers {
 
   int header_flag;           // 0 = item, 2 = xyz
   int flush_flag;            // 0 if no flush, 1 if flush every dump
+  int sort_flag;             // 1 if write in sorted order, 0 if not
 
   char *format_default;      // default format string
   char *format_user;         // format string set by user

src/dump_cfg.cpp

@@ -0,0 +1,403 @@
+/* ----------------------------------------------------------------------
+   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: Liang Wan (Chinese Academy of Sciences)
+------------------------------------------------------------------------- */
+
+#include "stdlib.h"
+#include "string.h"
+#include "dump_cfg.h"
+#include "atom.h"
+#include "domain.h"
+#include "comm.h"
+#include "modify.h"
+#include "compute.h"
+#include "input.h"
+#include "fix.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{INT,DOUBLE};  // same as in dump_custom.cpp
+
+/* ---------------------------------------------------------------------- */
+
+DumpCFG::DumpCFG(LAMMPS *lmp, int narg, char **arg) :
+  DumpCustom(lmp, narg, arg)
+{
+  if (narg < 10 ||
+      strcmp(arg[5],"id") != 0 || strcmp(arg[6],"type") != 0 ||
+      strcmp(arg[7],"xs") != 0 || strcmp(arg[8],"ys") != 0 ||
+      strcmp(arg[9],"zs") != 0)
+    error->all("Dump cfg arguments must start with 'id type xs ys zs'");
+
+  ntypes = atom->ntypes;
+  typenames = NULL;
+
+  // arrays for data rearrangement
+
+  recvcounts = new int[nprocs];
+  displs = new int[nprocs];
+  tags = NULL;
+  rbuf = NULL;
+  nchosen = nlines = 0;
+
+  // setup auxiliary property name strings
+  // convert 'X_ID[m]' (X=c,f,v) to 'ID_m'
+
+  if (narg > 10) auxname = new char*[narg-10];
+  else auxname = NULL;
+
+  int i = 0;
+  for (int iarg = 10; iarg < narg; iarg++, i++) {
+    if (strncmp(arg[iarg],"c_",2) == 0 ||
+	strncmp(arg[iarg],"f_",2) == 0 ||
+	strncmp(arg[iarg],"v_",2) == 0) {
+      int n = strlen(arg[iarg]);
+      char *suffix = new char[n];
+      strcpy(suffix,&arg[iarg][2]);
+
+      char *ptr = strchr(suffix,'[');
+      if (ptr) {
+	if (suffix[strlen(suffix)-1] != ']')
+	  error->all("Invalid keyword in dump cfg command");
+	*ptr = '\0';
+	*(ptr+2) = '\0';
+	auxname[i] = new char[strlen(suffix) + 3];
+	strcpy(auxname[i],suffix);
+	strcat(auxname[i],"_");
+	strcat(auxname[i],ptr+1);
+      } else {
+	auxname[i] = new char[strlen(suffix) + 1];
+	strcpy(auxname[i],suffix);
+      }
+
+      delete [] suffix;
+
+    } else {
+      auxname[i] = new char[strlen(arg[iarg]) + 1];
+      strcpy(auxname[i],arg[iarg]);
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+DumpCFG::~DumpCFG()
+{
+  if (typenames) {
+    for (int i = 1; i <= ntypes; i++) delete [] typenames[i];
+    delete [] typenames;
+  }
+
+  delete [] recvcounts;
+  delete [] displs;
+
+  if (tags) memory->sfree(tags);
+  if (rbuf) memory->destroy_2d_double_array(rbuf);
+
+  if (auxname) {
+    for (int i = 0; i < nfield-5; i++) delete [] auxname[i];
+    delete [] auxname;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCFG::init()
+{
+  if (multifile == 0)
+    error->all("Dump in CFG format requires one snapshot per file");
+
+  if (typenames == NULL) {
+    typenames = new char*[ntypes+1];
+    for (int itype = 1; itype <= ntypes; itype++) {
+      typenames[itype] = new char[3];
+      strcpy(typenames[itype],"C");
+    }
+    if (comm->me == 0)
+      error->warning("All element names have been set to 'C' for dump cfg");
+  }
+
+  // setup format strings
+
+  delete [] format;
+  char *str;
+  if (format_user) str = format_user;
+  else str = format_default;
+
+  int n = strlen(str) + 1;
+  format = new char[n];
+  strcpy(format,str);
+
+  // tokenize the format string and add space at end of each format element
+
+  char *ptr;
+  for (int i = 0; i < size_one; i++) {
+    if (i == 0) ptr = strtok(format," \0");
+    else ptr = strtok(NULL," \0");
+    if (ptr == NULL) error->all("Dump cfg user format string error");
+    delete [] vformat[i];
+    vformat[i] = new char[strlen(ptr) + 2];
+    strcpy(vformat[i],ptr);
+    vformat[i] = strcat(vformat[i]," ");
+  }
+  if (strtok(NULL," \0"))
+    error->all("Dump cfg user format string error");
+
+  // find current ptr for each compute,fix,variable
+  // check that fix frequency is acceptable
+
+  int icompute;
+  for (int i = 0; i < ncompute; i++) {
+    icompute = modify->find_compute(id_compute[i]);
+    if (icompute < 0) error->all("Could not find dump cfg compute ID");
+    compute[i] = modify->compute[icompute];
+  }
+
+  int ifix;
+  for (int i = 0; i < nfix; i++) {
+    ifix = modify->find_fix(id_fix[i]);
+    if (ifix < 0) error->all("Could not find dump cfg fix ID");
+    fix[i] = modify->fix[ifix];
+    if (nevery % modify->fix[ifix]->peratom_freq)
+      error->all("Dump cfg and fix not computed at compatible times");
+  }
+
+  int ivariable;
+  for (int i = 0; i < nvariable; i++) {
+    ivariable = input->variable->find(id_variable[i]);
+    if (ivariable < 0) error->all("Could not find dump cfg variable name");
+    variable[i] = ivariable;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCFG::write_header(int n)
+{
+  if (me == 0 || multiproc) {
+    fprintf(fp,"Number of particles = %d\n", n);
+    fprintf(fp,"A = 1.0 Angstrom (basic length-scale)\n");
+    fprintf(fp,"H0(1,1) = %g A\n",domain->xprd);
+    fprintf(fp,"H0(1,2) = 0 A \n");
+    fprintf(fp,"H0(1,3) = 0 A \n");
+    fprintf(fp,"H0(2,1) = %g A \n",domain->xy);
+    fprintf(fp,"H0(2,2) = %g A\n",domain->yprd);
+    fprintf(fp,"H0(2,3) = 0 A \n");
+    fprintf(fp,"H0(3,1) = %g A \n",domain->xz);
+    fprintf(fp,"H0(3,2) = %g A \n",domain->yz);
+    fprintf(fp,"H0(3,3) = %g A\n",domain->zprd);
+    fprintf(fp,".NO_VELOCITY.\n");
+    fprintf(fp,"entry_count = %d\n",nfield-2);
+    for (int i = 0; i < nfield-5; i++)
+      fprintf(fp,"auxiliary[%d] = %s\n",i,auxname[i]);
+  }
+
+  // calculate total # of data lines to be written on a writing proc
+
+  if (multiproc) nchosen = nmine;
+  else MPI_Reduce(&nmine,&nchosen,1,MPI_INT,MPI_SUM,0,world);
+
+  // allocate memory needed for data rearrangement on writing proc(s)
+
+  if (multiproc || me == 0) {
+    if (rbuf) memory->destroy_2d_double_array(rbuf);
+    rbuf = memory->create_2d_double_array(nchosen,size_one,"dump:rbuf");
+  }
+
+  // create a sorted list of atom IDs on writing proc(s) if necessary
+
+  if (sort_flag) create_sorted_tags();
+}
+
+/* ----------------------------------------------------------------------
+   write data lines to file in a block-by-block style
+   write head of block (mass & element name) only if has atoms of the type
+------------------------------------------------------------------------- */
+
+void DumpCFG::write_data(int n, double *buf)
+{
+  int i,j,m,itype;
+  int tag_i,index;
+  double *mass = atom->mass;
+
+  // if sort flag is off, transfer data in buf to rbuf directly
+  // if write by proc 0, transfer chunk by chunk
+
+  if (sort_flag == 0) {
+    for (i = 0, m = 0; i < n; i++) {
+      for (j = 0; j < size_one; j++)
+	rbuf[nlines][j] = buf[m++];
+      nlines++;
+    }
+
+    //  write data lines in rbuf to file after transfer is done
+
+    if (nlines == nchosen) {
+      for (itype = 1; itype <= ntypes; itype++) {
+	for (i = 0; i < nchosen; i++)
+	  if (rbuf[i][1] == itype) break;
+	if (i < nchosen) {
+	  fprintf(fp,"%g\n",mass[itype]);
+	  fprintf(fp,"%s\n",typenames[itype]);
+	  for (; i < nchosen; i++) {
+	    if (rbuf[i][1] == itype) {
+	      for (j = 2; j < size_one; j++) {
+		if (vtype[j] == INT)
+		  fprintf(fp,vformat[j],static_cast<int> (rbuf[i][j]));
+		else fprintf(fp,vformat[j],rbuf[i][j]);
+	      }
+	      fprintf(fp,"\n");
+	    }
+	  }
+	}
+      }
+      nlines = 0;
+    }
+
+  // if sort flag is on, transfer data lines in buf to rbuf &
+  //   rearrange them in a sorted order
+  // if write by proc 0, transfer chunk by chunk
+
+  } else {
+
+    // sort data lines:
+    // find the index of the atom ID of each data line within the
+    //   sorted atom IDs array
+    // transfer the data line to rbuf with its location according
+    //   to this index
+
+    for (i = 0, m = 0; i < n; i++) {
+      tag_i = static_cast<int>(buf[m]);
+      for (j = 0; j < nchosen; j++) {
+	if (tag_i == tags[j]) {
+	  index = j;
+	  break;
+	}
+      }
+      for (j = 0; j < size_one; j++)
+	rbuf[index][j] = buf[m++];
+    }
+
+    //  write data lines in rbuf to file after transfer is done
+
+    nlines += n;
+    if (nlines == nchosen) {
+      for (itype = 1; itype <= ntypes; itype++) {
+	for (i = 0; i < nchosen; i++)
+	  if (rbuf[i][1] == itype) break;
+	if (i < nchosen) {
+	  fprintf(fp,"%g\n",mass[itype]);
+	  fprintf(fp,"%s\n",typenames[itype]);
+	  for (; i < nchosen; i++) {
+	    if (rbuf[i][1] == itype) {
+	      for (j = 2; j < size_one; j++) {
+		if (vtype[j] == INT)
+		  fprintf(fp,vformat[j],static_cast<int>(rbuf[i][j]));
+		else fprintf(fp,vformat[j],rbuf[i][j]);
+	      }
+	      fprintf(fp,"\n");
+	    }
+	  }
+	}
+      }
+      nlines = 0;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   create a sorted list of atom IDs on writing proc(s)
+------------------------------------------------------------------------- */
+
+void DumpCFG::create_sorted_tags()
+{
+  int index;
+  int *mytags;
+  int *tag = atom->tag;
+  int nlocal = atom->nlocal;
+
+  index = 0;
+
+  // if write by multiproc, each proc has its own atom IDs collected
+
+  if (multiproc) {
+    if (tags) memory->sfree(tags);
+    tags = (int *) memory->smalloc(nchosen*sizeof(int),"dump:tags");
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) tags[index++] = tag[i];
+
+  // if write by proc 0, gather atom IDs of all the chosen atoms on proc 0
+
+  } else {
+    mytags = (int *) memory->smalloc(nmine*sizeof(int),"dump:mytags");
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) mytags[index++] = tag[i];
+    MPI_Gather(&nmine,1,MPI_INT,recvcounts,1,MPI_INT,0,world);
+    if (me == 0) {
+      if (tags) memory->sfree(tags);
+      tags = (int *) memory->smalloc(nchosen*sizeof(int),"dump:tags");
+      displs[0] = 0;
+      for (int iproc = 1; iproc < nprocs; iproc++)
+	displs[iproc] = displs[iproc-1] + recvcounts[iproc-1];
+    }
+    MPI_Gatherv(mytags,nmine,MPI_INT,tags,recvcounts,displs,MPI_INT,0,world);
+    memory->sfree(mytags);
+  }
+
+  if (multiproc || me == 0) qsort(tags,nchosen,sizeof(int),tag_compare);
+}
+
+/* ----------------------------------------------------------------------
+   compare function for qsort
+------------------------------------------------------------------------- */
+
+int DumpCFG::tag_compare(const void *itag, const void *jtag)
+{
+  int tag_i = *((int *) itag);
+  int tag_j = *((int *) jtag);
+
+  if (tag_i < tag_j) return -1;
+  else if (tag_i > tag_j) return 1;
+  return 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int DumpCFG::modify_param2(int narg, char **arg)
+{
+  if (strcmp(arg[0],"element") == 0) {
+    if (narg != ntypes+1)
+      error->all("Dump modify element names do not match atom types");
+
+    if (typenames) {
+      for (int i = 1; i <= ntypes; i++) delete [] typenames[i];
+      delete [] typenames;
+      typenames = NULL;
+    }
+
+    typenames = new char*[ntypes+1];
+    for (int itype = 1; itype <= ntypes; itype++) {
+      typenames[itype] = new char[3];
+      if (strlen(arg[itype]) >= 3)
+	error->all("Illegal chemical element names");
+      strcpy(typenames[itype],arg[itype]);
+    }
+    return ntypes+1;
+
+  } else return 0;
+}

src/dump_cfg.h

@@ -0,0 +1,47 @@
+/* ----------------------------------------------------------------------
+   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 DUMP_CFG_H
+#define DUMP_CFG_H
+
+#include "dump_custom.h"
+
+namespace LAMMPS_NS {
+
+class DumpCFG : public DumpCustom {
+ public:
+  DumpCFG(class LAMMPS *, int, char **);
+  ~DumpCFG();
+  void init();
+  void write_header(int);
+  void write_data(int, double *);
+  int modify_param2(int, char **);
+
+ private:
+  int ntypes;                // # of atom types
+  char **typenames;	     // array of element names for each type
+  char **auxname;            // name strings of auxiliary properties
+  int nchosen;               // # of lines to be written on a writing proc
+  int nlines;                // # of lines transferred from buf to rbuf
+  int *recvcounts;           // # of data (tag info) received from all procs
+  int *displs;               // displacement of data (tag info) within buffer
+  int *tags;                 // atom IDs of all the chosen atoms collected
+  double **rbuf;             // buf of data lines for data lines rearrangement
+
+  void create_sorted_tags();
+  static int tag_compare(const void *, const void *);
+};
+
+}
+
+#endif

src/dump_custom.cpp

@@ -1365,6 +1365,12 @@ int DumpCustom::modify_param(int narg, char **arg)
 
     nthresh++;
     return 4;
+
+  // pass along params to child class
+
+  } else {
+    int n = modify_param2(narg,arg);
+    return n;
   }
 
   return 0;

src/dump_custom.h

@@ -21,11 +21,11 @@ namespace LAMMPS_NS {
 class DumpCustom : public Dump {
  public:
   DumpCustom(class LAMMPS *, int, char **);
-  ~DumpCustom();
+  virtual ~DumpCustom();
-  void init();
+  virtual void init();
   double memory_usage();
 
- private:
+ protected:
   int nevery;                // dump frequency to check Fix against
   int iregion;               // -1 if no region, else which region
   int nthresh;               // # of defined threshholds
@@ -64,16 +64,17 @@ class DumpCustom : public Dump {
 
   // private methods
 
-  void write_header(int);
+  virtual void write_header(int);
   int count();
   int pack();
-  void write_data(int, double *);
+  virtual void write_data(int, double *);
 
   void parse_fields(int, char **);
   int add_compute(char *);
   int add_fix(char *);
   int add_variable(char *);
   int modify_param(int, char **);
+  virtual int modify_param2(int, char **) {return 0;}
 
   typedef void (DumpCustom::*FnPtrHeader)(int);
   FnPtrHeader header_choice;           // ptr to write header functions

src/style.h

@@ -135,6 +135,7 @@ ComputeStyle(temp/sphere,ComputeTempSphere)
 
 #ifdef DumpInclude
 #include "dump_atom.h"
+#include "dump_cfg.h"
 #include "dump_custom.h"
 #include "dump_dcd.h"
 #include "dump_xyz.h"
@@ -142,6 +143,7 @@ ComputeStyle(temp/sphere,ComputeTempSphere)
 
 #ifdef DumpClass
 DumpStyle(atom,DumpAtom)
+DumpStyle(cfg,DumpCFG)
 DumpStyle(custom,DumpCustom)
 DumpStyle(dcd,DumpDCD)
 DumpStyle(xyz,DumpXYZ)