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

src/fix_property_atom.cpp

@@ -0,0 +1,348 @@
+/* ----------------------------------------------------------------------
+   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 "stdlib.h"
+#include "string.h"
+#include "fix_property_atom.h"
+#include "atom.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+enum{MOLECULE,INTEGER,DOUBLE};
+
+  // NOTE: how to prevent atom style being changed after this fix defined?
+  //       e.g. what if replicate happens, will it wipe out molecule_flag
+  //       will it copy ivectors in atom.h ??
+  //       maybe setting of atom->mol_flag should happen in fix init
+
+/* ---------------------------------------------------------------------- */
+
+FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 4) error->all(FLERR,"Illegal fix property/atom command");
+
+  restart_peratom = 1;
+
+  int iarg = 3;
+  nvalue = narg-iarg;
+  style = new int[nvalue];
+
+  molecule_flag = 0;
+
+  nvalue = 0;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"mol") == 0) {
+      if (atom->molecule_flag)
+        error->all(FLERR,"Fix property/atom mol when atom_style "
+                   "already has molecule attribute");
+      if (molecule_flag) 
+        error->all(FLERR,"Fix property/atom cannot specify mol twice");
+      style[nvalue] = MOLECULE;
+      atom->molecule_flag = molecule_flag = 1;
+      nvalue++;
+    } else if (strstr(arg[iarg],"i_") == arg[iarg]) {
+      style[nvalue] = INTEGER;
+      nvalue++;
+    } else if (strstr(arg[iarg],"d_") == arg[iarg]) {
+      style[nvalue] = DOUBLE;
+      nvalue++;
+    } else break;
+
+    iarg++;
+  }
+
+  // optional args
+
+  border = 0;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"border") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix property/atom command");
+      if (strcmp(arg[iarg+1],"no") == 0) border = 0;
+      else if (strcmp(arg[iarg+1],"yes") == 0) border = 1;
+      else error->all(FLERR,"Illegal fix property/atom command");
+    } else error->all(FLERR,"Illegal fix property/atom command");
+  }
+
+  // perform initial allocation of atom-based array
+  // register with Atom class
+
+  nmax_old = 0;
+  grow_arrays(atom->nmax);
+  atom->add_callback(0);
+  atom->add_callback(1);
+  if (border) atom->add_callback(2);
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixPropertyAtom::~FixPropertyAtom()
+{
+  // unregister callbacks to this fix from Atom class
+
+  atom->delete_callback(id,0);
+  atom->delete_callback(id,1);
+  if (border) atom->delete_callback(id,2);
+
+  // deallocate per-atom vectors in Atom class
+  // set ptrs to NULL, so they no longer exist for Atom class
+
+  for (int m = 0; m < nvalue; m++) {
+    if (style[m] == MOLECULE) {
+      atom->molecule_flag = 0;
+      memory->destroy(atom->molecule);
+      atom->molecule = NULL;
+    } else if (style[m] == INTEGER) {
+    } else if (style[m] == DOUBLE) {
+    }
+  }
+
+  delete [] style;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixPropertyAtom::setmask()
+{
+  int mask = 0;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPropertyAtom::init()
+{
+  // check again if atom style already defines molecule vector
+  //   based on molecular setting
+  // this could happen if script changed atom_style
+  //   after this fix was already in place
+
+  if (molecule_flag && atom->molecular)
+    error->all(FLERR,"Fix property/atom mol when atom_style "
+               "already has molecule attribute");
+}
+
+/* ----------------------------------------------------------------------
+   unpack section of data file
+------------------------------------------------------------------------- */
+
+void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf)
+{
+  int j,m,tagdata;
+  char *next;
+
+  next = strchr(buf,'\n');
+  *next = '\0';
+  int nwords = atom->count_words(buf);
+  *next = '\n';
+
+  if (nwords != nvalue+1) {
+    char str[128];
+    sprintf(str,"Incorrect %s format in data file",keyword);
+    error->all(FLERR,str);
+  }
+
+  char **values = new char*[nwords];
+
+  // loop over lines of atom velocities
+  // tokenize the line into values
+  // if I own atom tag, unpack its values
+
+  int map_tag_max = atom->map_tag_max;
+
+  for (int i = 0; i < n; i++) {
+    next = strchr(buf,'\n');
+
+    values[0] = strtok(buf," \t\n\r\f");
+    for (j = 1; j < nwords; j++)
+      values[j] = strtok(NULL," \t\n\r\f");
+
+    tagdata = atoi(values[0]);
+    if (tagdata <= 0 || tagdata > map_tag_max) {
+      char str[128];
+      sprintf(str,"Invalid atom ID in %s section of data file",keyword);
+      error->one(FLERR,str);
+    }
+
+    // assign words in line to per-atom vectors
+
+    if ((m = atom->map(tagdata)) >= 0) {
+      for (j = 0; j < nvalue; j++) {
+        if (style[j] == MOLECULE) atom->molecule[m] = atoi(values[j+1]);
+        else if (style[j] == INTEGER) atom->molecule[m] = atoi(values[j+1]);
+        else if (style[j] == DOUBLE) atom->molecule[m] = atof(values[j+1]);
+      }
+    }
+
+    buf = next + 1;
+  }
+
+  delete [] values;
+}
+
+/* ---------------------------------------------------------------------- */
+
+bigint FixPropertyAtom::read_data_skip_lines(char *keyword)
+{
+  return atom->natoms;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double FixPropertyAtom::memory_usage()
+{
+  double bytes = 0.0;
+  for (int m = 0; m < nvalue; m++) {
+    if (style[m] == MOLECULE) bytes = atom->nmax * sizeof(int);
+    else if (style[m] == INTEGER) bytes = atom->nmax * sizeof(int);
+    else if (style[m] == DOUBLE) bytes = atom->nmax * sizeof(double);
+  }
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+   allocate atom-based arrays
+   initialize new values to 0,
+   since AtomVec class won't do it as atoms are added, 
+   e.g. in create_atom() or data_atom()
+------------------------------------------------------------------------- */
+
+void FixPropertyAtom::grow_arrays(int nmax)
+{
+  for (int m = 0; m < nvalue; m++) {
+    if (style[m] == MOLECULE) {
+      memory->grow(atom->molecule,nmax,"atom:molecule");
+      size_t nbytes = (nmax-nmax_old) * sizeof(int);
+      memset(&atom->molecule[nmax_old],0,nbytes);
+    } else if (style[m] == INTEGER) {
+      memory->grow(atom->molecule,nmax,"atom:molecule");
+      size_t nbytes = (nmax-nmax_old) * sizeof(int);
+      memset(&atom->molecule[nmax_old],0,nbytes);
+    } else if (style[m] == DOUBLE) {
+      memory->grow(atom->molecule,nmax,"atom:molecule");
+      size_t nbytes = (nmax-nmax_old) * sizeof(double);
+      memset(&atom->molecule[nmax_old],0,nbytes);
+    }
+  }
+
+  nmax_old = nmax;
+}
+
+/* ----------------------------------------------------------------------
+   copy values within local atom-based arrays
+------------------------------------------------------------------------- */
+
+void FixPropertyAtom::copy_arrays(int i, int j, int delflag)
+{
+  for (int m = 0; m < nvalue; m++) {
+    if (style[m] == MOLECULE)
+      atom->molecule[j] = atom->molecule[i];
+    else if (style[m] == INTEGER)
+      atom->molecule[j] = atom->molecule[i];
+    else if (style[m] == DOUBLE)
+      atom->molecule[j] = atom->molecule[i];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based array for exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixPropertyAtom::pack_exchange(int i, double *buf)
+{
+  for (int m = 0; m < nvalue; m++) {
+    if (style[m] == MOLECULE) buf[m] = atom->molecule[i];
+    else if (style[m] == INTEGER) buf[m] = atom->molecule[i];
+    else if (style[m] == DOUBLE) buf[m] = atom->molecule[i];
+  }
+  return nvalue;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values in local atom-based array from exchange with another proc
+------------------------------------------------------------------------- */
+
+int FixPropertyAtom::unpack_exchange(int nlocal, double *buf)
+{
+  for (int m = 0; m < nvalue; m++) {
+    if (style[m] == MOLECULE) 
+      atom->molecule[nlocal] = static_cast<int> (buf[m]);
+    else if (style[m] == INTEGER) 
+      atom->molecule[nlocal] = static_cast<int> (buf[m]);
+    else if (style[m] == DOUBLE)
+      atom->molecule[nlocal] = buf[m];
+  }
+  return nvalue;
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based arrays for restart file
+------------------------------------------------------------------------- */
+
+int FixPropertyAtom::pack_restart(int i, double *buf)
+{
+  buf[0] = nvalue+1;
+  for (int m = 1; m <= nvalue; m++) {
+    if (style[m] == MOLECULE) buf[m] = atom->molecule[i];
+    else if (style[m] == INTEGER) buf[m] = atom->molecule[i];
+    else if (style[m] == DOUBLE) buf[m] = atom->molecule[i];
+  }
+  return nvalue+1;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values from atom->extra array to restart the fix
+------------------------------------------------------------------------- */
+
+void FixPropertyAtom::unpack_restart(int nlocal, int nth)
+{
+  double **extra = atom->extra;
+
+  // skip to Nth set of extra values
+
+  int m = 0;
+  for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
+  m++;
+
+  for (int i = 0; i < nvalue; i++) {
+    if (style[i] == MOLECULE) 
+      atom->molecule[nlocal] = static_cast<int> (extra[nlocal][m++]);
+    else if (style[m] == INTEGER) 
+      atom->molecule[nlocal] = static_cast<int> (extra[nlocal][m++]);
+    else if (style[m] == DOUBLE)
+      atom->molecule[nlocal] = extra[nlocal][m++];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   maxsize of any atom's restart data
+------------------------------------------------------------------------- */
+
+int FixPropertyAtom::maxsize_restart()
+{
+  return nvalue+1;
+}
+
+/* ----------------------------------------------------------------------
+   size of atom nlocal's restart data
+------------------------------------------------------------------------- */
+
+int FixPropertyAtom::size_restart(int nlocal)
+{
+  return nvalue+1;
+}

src/fix_property_atom.h

@@ -0,0 +1,60 @@
+/* ----------------------------------------------------------------------
+   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(property/atom,FixPropertyAtom)
+
+#else
+
+#ifndef LMP_FIX_PROPERTY_ATOM_H
+#define LMP_FIX_PROPERTY_ATOM_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixPropertyAtom : public Fix {
+ public:
+  FixPropertyAtom(class LAMMPS *, int, char **);
+  ~FixPropertyAtom();
+  int setmask();
+  void init();
+
+  void read_data_section(char *, int, char *);
+  bigint read_data_skip_lines(char *);
+
+  void grow_arrays(int);
+  void copy_arrays(int, int, int);
+  int pack_exchange(int, double *);
+  int unpack_exchange(int, double *);
+  int pack_restart(int, double *);
+  void unpack_restart(int, int);
+  int size_restart(int);
+  int maxsize_restart();
+  double memory_usage();
+
+ private:
+  int nvalue,border,molecule_flag;
+  int *style;
+  int nmax_old;         // length of peratom arrays the last time they grew
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+*/

src/modify.cpp

@@ -634,15 +634,16 @@ int Modify::min_reset_ref()
 
 void Modify::add_fix(int narg, char **arg, char *suffix)
 {
-  const char *exceptions[NEXCEPT] = {"GPU","OMP","atom/property","cmap"};
-
   if (narg < 3) error->all(FLERR,"Illegal fix command");
 
   // cannot define fix before box exists unless style is in exception list
   // don't like this way of checking for exceptions by adding fixes to list,
   //   but can't think of better way
   // too late if instantiate fix, then check flag set in fix constructor,
-  // since some fixes access domain settings in their constructor
+  //   since some fixes access domain settings in their constructor
+  // change NEXCEPT above when add new fix to this list
+
+  const char *exceptions[NEXCEPT] = {"GPU","OMP","property/atom","cmap"};
 
   if (domain->box_exist == 0) {
     int m;

src/read_data.cpp

@@ -124,10 +124,13 @@ void ReadData::command(int narg, char **arg)
       fix_index[nfix] = modify->find_fix(arg[iarg+1]);
       if (fix_index[nfix] < 0)
         error->all(FLERR,"Fix ID for read_data does not exist");
-      int n = strlen(arg[iarg+2]) + 1;
+      if (strcmp(arg[iarg+2],"NULL") == 0) fix_header[nfix] = NULL;
-      fix_header[nfix] = new char[n];
+      else {
-      strcpy(fix_header[nfix],arg[iarg+2]);
+        int n = strlen(arg[iarg+2]) + 1;
-      n = strlen(arg[iarg+3]) + 1;
+        fix_header[nfix] = new char[n];
+        strcpy(fix_header[nfix],arg[iarg+2]);
+      }
+      int n = strlen(arg[iarg+3]) + 1;
       fix_section[nfix] = new char[n];
       strcpy(fix_section[nfix],arg[iarg+3]);
       nfix++;
@@ -198,10 +201,8 @@ void ReadData::command(int narg, char **arg)
 
     if (nfix) {
       for (n = 0; n < nfix; n++)
-        if (strstr(line,fix_section[n])) {
+        if (strcmp(keyword,fix_section[n]) == 0) {
-          bigint nlines =
+          fix(n,keyword);
-            modify->fix[fix_index[n]]->read_data_skip_lines(keyword);
-          fix(n,keyword,nlines);
           parse_keyword(0,1);
           break;
         }
@@ -445,11 +446,13 @@ void ReadData::header(int flag)
     // if fix matches, continue to next header line
 
     if (nfix) {
-      for (n = 0; n < nfix; n++)
+      for (n = 0; n < nfix; n++) {
+        if (!fix_header[n]) continue;
         if (strstr(line,fix_header[n])) {
           modify->fix[fix_index[n]]->read_data_header(line);
           break;
         }
+      }
       if (n < nfix) continue;
     }
 
@@ -1089,16 +1092,18 @@ void ReadData::impropercoeffs(int which)
    n = index of fix
 ------------------------------------------------------------------------- */
 
-void ReadData::fix(int ifix, char *line, bigint nlines)
+void ReadData::fix(int ifix, char *keyword)
 {
   int nchunk,eof;
 
+  bigint nlines = modify->fix[ifix]->read_data_skip_lines(keyword);
+
   bigint nread = 0;
   while (nread < nlines) {
     nchunk = MIN(nlines-nread,CHUNK);
     eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
     if (eof) error->all(FLERR,"Unexpected end of data file");
-    modify->fix[ifix]->read_data_section(line,nchunk,buffer);
+    modify->fix[ifix]->read_data_section(keyword,nchunk,buffer);
     nread += nchunk;
   }
 }
@@ -1141,10 +1146,8 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
 
     if (nfix) {
       for (i = 0; i < nfix; i++) {
-        printf("LINE SECTION %s %s\n",line,fix_section[i]);
+        if (strcmp(keyword,fix_section[i]) == 0) {
-        if (strstr(line,fix_section[i])) {
           int n = modify->fix[fix_index[i]]->read_data_skip_lines(keyword);
-          printf("NLINES SKIP %d\n",n);
           skip_lines(n);
           parse_keyword(0,0);
           break;
@@ -1512,6 +1515,8 @@ void ReadData::parse_keyword(int first, int flag)
 
 /* ----------------------------------------------------------------------
    proc 0 reads N lines from file
+   NOTE: needs to be called with bigint in some cases
+         if called with int, will it be promoted to bigint?
 ------------------------------------------------------------------------- */
 
 void ReadData::skip_lines(int n)

src/read_data.h

@@ -78,7 +78,7 @@ class ReadData : protected Pointers {
   void dihedralcoeffs(int);
   void impropercoeffs(int);
 
-  void fix(int, char *, bigint);
+  void fix(int, char *);
 };
 
 }

src/replicate.cpp

@@ -82,7 +82,7 @@ void Replicate::command(int narg, char **arg)
   // maxmol = largest molecule tag across all existing atoms
 
   int maxmol = 0;
-  if (atom->molecular) {
+  if (atom->molecule_flag) {
     for (i = 0; i < atom->nlocal; i++) maxmol = MAX(atom->molecule[i],maxmol);
     int maxmol_all;
     MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_INT,MPI_MAX,world);
@@ -127,7 +127,7 @@ void Replicate::command(int narg, char **arg)
   // if molecular and N > MAXTAGINT, error
   // if atomic and new N > MAXTAGINT, turn off tags for existing and new atoms
   // new system cannot exceed MAXBIGINT
-  // change these 2 to MAXTAGINT when allow tagint = bigint
+  // NOTE: change these 2 to MAXTAGINT when allow tagint = bigint
 
   if (atom->molecular && 
       (nrep*old->natoms < 0 || nrep*old->natoms > MAXSMALLINT))
@@ -315,7 +315,7 @@ void Replicate::command(int narg, char **arg)
               atom->tag[i] += atom_offset;
               atom->image[i] = image;
 
-              if (atom->molecular) {
+              if (atom->molecule_flag) {
                 if (atom->molecule[i] > 0)
                   atom->molecule[i] += mol_offset;
                 if (atom->avec->bonds_allow)