// ATC headers
#include "ATC_Method.h"
#include "LammpsInterface.h"
#include "FE_Engine.h"
#include "Array.h"
#include "Array2D.h"
#include "ATC_Error.h"
#include "Function.h"
#include "PrescribedDataManager.h"
#include "TimeIntegrator.h"
#include "PhysicsModel.h"
#include "PerAtomQuantityLibrary.h"
#include "TransferLibrary.h"
#include "KernelFunction.h"
#include "Utility.h"
#include "FieldManager.h"

#include <fstream>
#include <sstream>
#include <iostream>

using ATC_Utility::sgn;
using ATC_Utility::to_string;
using ATC_Utility::is_dbl;

using std::stringstream;
using std::ifstream;
using std::ofstream;
using std::string;
using std::map;
using std::set;
using std::vector;
using std::pair;

namespace ATC {

  ATC_Method::ATC_Method(string groupName, double ** & perAtomArray, LAMMPS_NS::Fix * thisFix) :
    nodalAtomicVolume_(nullptr),
    needReset_(true),
    lammpsInterface_(LammpsInterface::instance()),
    interscaleManager_(this),
    timeFilterManager_(this),
    integrateInternalAtoms_(false),
    atomTimeIntegrator_(nullptr),
    ghostManager_(this),
    feEngine_(nullptr),
    initialized_(false),
    meshDataInitialized_(false),
    localStep_(0),
    sizeComm_(8), // 3 positions + 1 material id * 2 for output
    atomCoarseGrainingPositions_(nullptr),
    atomGhostCoarseGrainingPositions_(nullptr),
    atomProcGhostCoarseGrainingPositions_(nullptr),
    atomReferencePositions_(nullptr),
    nNodes_(0),
    nsd_(lammpsInterface_->dimension()),
    xref_(nullptr),
    readXref_(false),
    needXrefProcessorGhosts_(false),
    trackDisplacement_(false),
    needsAtomToElementMap_(true),
    atomElement_(nullptr),
    atomGhostElement_(nullptr),
    internalElementSet_(""),
    atomMasses_(nullptr),
    atomPositions_(nullptr),
    atomVelocities_(nullptr),
    atomForces_(nullptr),
    parallelConsistency_(true),
    outputNow_(false),
    outputTime_(true),
    outputFrequency_(0),
    sampleFrequency_(0),
    sampleCounter_(0),
    peScale_(1./(lammpsInterface_->mvv2e())),
    keScale_(1.),
    scalarFlag_(0),
    vectorFlag_(0),
    sizeVector_(0),
    scalarVectorFreq_(0),
    sizePerAtomCols_(4),
    perAtomOutput_(nullptr),
    perAtomArray_(perAtomArray),
    extScalar_(0),
    extVector_(0),
    extList_(nullptr),
    thermoEnergyFlag_(0),
    atomVolume_(nullptr),
    atomicWeightsWriteFlag_(false),
    atomicWeightsWriteFrequency_(0),
    atomWeightType_(LATTICE),
    domainDecomposition_(REPLICATED_MEMORY),
    groupbit_(0),
    groupbitGhost_(0),
    needProcGhost_(false),
    groupTag_(groupName),
    nLocal_(0),
    nLocalTotal_(0),
    nLocalGhost_(0),
    atomToElementMapType_(LAGRANGIAN),
    atomToElementMapFrequency_(0),
    regionID_(-1),
    mdMassNormalization_(false),
    kernelBased_(false),
    kernelOnTheFly_(false),
    kernelFunction_(nullptr),
    bondOnTheFly_(false),
    accumulant_(nullptr),
    accumulantMol_(nullptr),
    accumulantMolGrad_(nullptr),
    accumulantWeights_(nullptr),
    accumulantInverseVolumes_(&invNodeVolumes_),
    accumulantBandwidth_(0),
    useRestart_(false),
    hasRefPE_(false),
    setRefPE_(false),
    setRefPEvalue_(false),
    refPEvalue_(0.),
    readRefPE_(false),
    nodalRefPotentialEnergy_(nullptr),
    simTime_(0.0),
    stepCounter_(0)
  {
    lammpsInterface_->print_msg_once("version "+version());
    lammpsInterface_->set_fix_pointer(thisFix);
    interscaleManager_.set_lammps_data_prefix();
    grow_arrays(lammpsInterface_->nmax());
    feEngine_ = new FE_Engine(lammpsInterface_->world());


    lammpsInterface_->create_compute_pe_peratom();
  }

  ATC_Method::~ATC_Method()
  {
    lammpsInterface_->destroy_2d_double_array(xref_);
    lammpsInterface_->destroy_2d_double_array(perAtomOutput_);
    if (atomTimeIntegrator_) delete atomTimeIntegrator_;
    if (feEngine_) delete feEngine_;
  }

  //--------------------------------------------------
  // pack_fields
  //   bundle all allocated field matrices into a list
  //   for output needs
  //--------------------------------------------------
  void ATC_Method::pack_fields(RESTART_LIST & data)
  {
    map<FieldName,int>::const_iterator field;
    for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
      FieldName thisField = field->first;
      string fieldName = field_to_string(thisField);
      string matrixName;
      // copy all fields from ATC_Method.h
      matrixName = "fields_" + fieldName;
      data[matrixName] = & fields_[thisField].set_quantity();
      matrixName = "dot_fields_" + fieldName;
      data[matrixName] = & dot_fields_[thisField].set_quantity();
      matrixName = "ddot_fields_" + fieldName;
      data[matrixName] = & ddot_fields_[thisField].set_quantity();
      matrixName = "dddot_fields_" + fieldName;
      data[matrixName] = & dddot_fields_[thisField].set_quantity();
      matrixName = "NodalAtomicFieldsRoc_" + fieldName;
      data[matrixName] = & nodalAtomicFieldsRoc_[thisField].set_quantity();
    }
  }

  //--------------------------------------------------
  // write_restart_file
  //   bundle matrices that need to be saved and call
  //   fe_engine to write the file
  //--------------------------------------------------
  void ATC_Method::write_restart_data(string fileName, RESTART_LIST & data)
  {
    pack_fields(data);
    feEngine_->write_restart_file(fileName,&data);
  }

  //--------------------------------------------------
  // read_restart_file
  //   bundle matrices that need to be saved and call
  //   fe_engine to write the file
  //--------------------------------------------------
  void ATC_Method::read_restart_data(string fileName, RESTART_LIST & data)
  {
    pack_fields(data);
    feEngine_->read_restart_file(fileName,&data);
  }

  //--------------------------------------------------
  // Interactions with LAMMPS fix commands
  // parse input command and pass on to finite element engine
  //   or physics specific transfers if necessary
  //   revert to physics-specific transfer if no command matches input
  // first keyword is unique to particular class
  // base class keyword matching must apply to ALL physics
  // order:  derived, base, owned objects
  //--------------------------------------------------
  bool ATC_Method::modify(int narg, char **arg)
  {
    int argIdx=0;
    bool match = false;

    // gateways to other modules e.g. extrinsic, control, mesh
    // pass off to fe engine
    if (strcmp(arg[argIdx],"mesh")==0) {
      match = feEngine_->modify(narg, arg);
      if (feEngine_->has_mesh()  && !meshDataInitialized_)
        this->initialize_mesh_data();
    }
    // pass off to time filter
    else if (strcmp(arg[argIdx],"filter")==0) {
      argIdx++;
      match = timeFilterManager_.modify(narg-argIdx,&arg[argIdx]);

        // consistentInitialization_ = false;
    }
    // pass off to kernel function manager
    else if (strcmp(arg[argIdx],"kernel")==0) {
      argIdx++;

      if (kernelFunction_) {
        //delete kernelFunction_;
        //resetKernelFunction_ = true;
      }
      kernelFunction_ = KernelFunctionMgr::instance()->function(&arg[argIdx],narg-argIdx);
      if (kernelFunction_) match = true;
      else ATC_Error("no matching kernel found");
      feEngine_->set_kernel(kernelFunction_);

      accumulantMol_=&kernelAccumulantMol_; // KKM add
      accumulantMolGrad_=&kernelAccumulantMolGrad_; // KKM add
    }
    // pass off to ghost manager
    else if (strcmp(arg[argIdx],"boundary_dynamics")==0) {
      argIdx++;
      match = ghostManager_.modify(narg-argIdx,&arg[argIdx]);
    }

    // parsing handled here
    else {
      if (strcmp(arg[argIdx],"parallel_consistency")==0) {
        argIdx++;
        //if (!kernelFunction_)          { throw ATC_Error("on_the_fly requires a kernel function"); }
        if (strcmp(arg[argIdx],"off")==0) parallelConsistency_ = false;
        else                              parallelConsistency_ = true;
        match = true;
      }
     /*! \page man_hardy_on_the_fly fix_modify AtC on_the_fly
        \section syntax
        fix_modify AtC on_the_fly <bond | kernel> <optional on | off> \n        - bond | kernel (keyword) = specifies on-the-fly calculation of bond or
kernel         matrix elements \n
        - on | off (keyword) =  activate or discontinue on-the-fly mode \n
        \section examples
        <TT> fix_modify AtC on_the_fly bond on </TT> \n        <TT> fix_modify AtC on_the_fly kernel </TT> \n
        <TT> fix_modify AtC on_the_fly kernel off </TT> \n
        \section description
        Overrides normal mode of pre-calculating and storing bond pair-to-node a
nd
        kernel atom-to-node matrices. If activated, will calculate elements of t
hese
        matrices during repeated calls of field computations (i.e. "on-the-fly") and not store them for
        future use.   \n        on flag is optional - if omitted, on_the_fly will be activated for the s
pecified
        matrix. Can be deactivated using off flag. \n
        \section restrictions
        Must be used with the hardy/field type of AtC fix
        ( see \ref man_fix_atc )
        \section related
        \section default
        By default, on-the-fly calculation is not active (i.e. off). However, code does a memory allocation check to determine if it can store all needed bond and kernel matrix ele ments. If this allocation fails, on-the-fly is activated. \n
      */

      else if (strcmp(arg[argIdx],"on_the_fly")==0) {
        argIdx++;
        //if (!kernelFunction_)          { throw ATC_Error("on_the_fly requires a kernel function"); }
        if (strcmp(arg[argIdx],"bond")==0) {
          argIdx++;
          bondOnTheFly_ = true;
          if (narg > argIdx && strcmp(arg[argIdx],"off")==0) bondOnTheFly_ = false;
        }
        else if (strcmp(arg[argIdx],"kernel")==0) {
          argIdx++;
          kernelOnTheFly_ = true;
          if (narg > argIdx && strcmp(arg[argIdx],"off")==0) kernelOnTheFly_ = false;
        }
        else { throw ATC_Error("unsupported on_the_fly type"); }
        match = true;
      }

      /*! \page man_output fix_modify AtC output
        \section syntax
        fix_modify AtC output <filename_prefix> <frequency>
        [text | full_text | binary | vector_components | tensor_components ]
        fix_modify AtC output index [step | time ]
        - filename_prefix (string) = prefix for data files
        - frequency (integer) = frequency of output in time-steps
        - options (keyword/s): \n
        text = creates text output of index, step and nodal variable values for unique nodes \n
        full_text = creates text output index, nodal id, step, nodal coordinates and nodal variable values for unique and image nodes \n
        binary = creates binary Ensight output \n
        vector_components = outputs vectors as scalar components \n
        tensor_components = outputs tensor as scalar components
        (use this for Paraview)\n

        \section examples
        <TT> fix_modify AtC output heatFE 100 </TT> \n
        <TT> fix_modify AtC output hardyFE 1 text tensor_components </TT> \n
        <TT> fix_modify AtC output hardyFE 10 text binary tensor_components </TT> \n
        <TT> fix_modify AtC output index step </TT> \n
        \section description
        Creates text and/or binary (Ensight, "gold" format) output of nodal/mesh data
        which is transfer/physics specific. Output indexed by step or time is possible.
        \section restrictions
        \section related
        see \ref man_fix_atc
        \section default
        no default format
        output indexed by time
      */
      else if (strcmp(arg[argIdx],"output")==0) {
        argIdx++;
      /*! \page man_output_nodeset fix_modify AtC output nodeset
        \section syntax
        fix_modify AtC output nodeset <nodeset_name> <operation>
        - nodeset_name (string) = name of nodeset to be operated on
        - operation (keyword/s): \n
        sum = creates nodal sum over nodes in specified nodeset \n
        \section examples
        <TT> fix_modify AtC output nodeset nset1 sum </TT> \n
        \section description
        Performs operation over the nodes belonging to specified nodeset
        and outputs resulting variable values to GLOBALS file.
        \section restrictions
        \section related
        see \ref man_fix_atc
        \section default
        none
      */
        if (strcmp(arg[argIdx],"nodeset")==0) {
          argIdx++;
          string nset = arg[argIdx++];
          if       (strcmp(arg[argIdx],"sum")==0) {
            argIdx++;
            string field = arg[argIdx];
            pair < string, FieldName >  id(nset,string_to_field(field));
            nsetData_[id] = NODESET_SUM;
            match = true;
          }
          else if (strcmp(arg[argIdx],"average")==0) {
            argIdx++;
            string field = arg[argIdx];
            pair < string, FieldName >  id(nset,string_to_field(field));
            nsetData_[id] = NODESET_AVERAGE;
            match = true;
          }
        }

      /*! \page man_boundary_integral fix_modify AtC output boundary_integral
        \section syntax
        fix_modify AtC output boundary_integral [field] faceset [name]
        - field (string) : name of hardy field
        - name (string)  : name of faceset
        \section examples
        <TT> fix_modify AtC output boundary_integral stress faceset loop1 </TT> \n
        \section description
        Calculates a surface integral of the given field dotted with the
        outward normal of the faces and puts output in the "GLOBALS" file
        \section restrictions
        Must be used with the hardy/field type of AtC fix
        ( see \ref man_fix_atc )
        \section related
        \section default
        none
      */

      /*! \page man_contour_integral fix_modify AtC output contour_integral
        \section syntax
        fix_modify AtC output contour_integral [field] faceset [name] <axis [x | y | z
]>
        - field (string) : name of hardy field
        - name (string)  : name of faceset
        - axis (string)  : x or y or z
        \section examples
        <TT> fix_modify AtC output contour_integral stress faceset loop1 </TT> \n
        \section description
        Calculates a surface integral of the given field dotted with the
        outward normal of the faces and puts output in the "GLOBALS" file
        \section restrictions
        Must be used with the hardy/field type of AtC fix
        ( see \ref man_fix_atc )
        \section related
        \section default
        none
      */

        else if ( (strcmp(arg[argIdx],"boundary_integral")==0)
               || (strcmp(arg[argIdx],"contour_integral")==0) ) {
          FacesetIntegralType type = BOUNDARY_INTEGRAL;
          if  (strcmp(arg[argIdx],"contour_integral")==0)
                              type = CONTOUR_INTEGRAL;
          argIdx++;
          string field(arg[argIdx++]);
          if(strcmp(arg[argIdx],"faceset")==0) {
            argIdx++;
            string name(arg[argIdx++]);
            pair <string,string> pair_name(name,field);
            fsetData_[pair_name] = type;
            match = true;
          }
        } // end "boundary_integral" || "contour_integral"

      /*! \page man_output_elementset fix_modify AtC output elementset
        \section syntax
        fix_modify AtC output volume_integral <eset_name> <field>
        - set_name (string) = name of elementset to be integrated over
        - fieldname (string) = name of field to integrate
        csum = creates nodal sum over nodes in specified nodeset \n
        \section examples
        <TT> fix_modify AtC output eset1 mass_density </TT> \n
        \section description
        Performs volume integration of specified field over elementset
        and outputs resulting variable values to GLOBALS file.
        \section restrictions
        \section related
        see \ref man_fix_atc
        \section default
        none
      */

        else if ( (strcmp(arg[argIdx],"volume_integral")==0) ) {
          argIdx++;
          string name(arg[argIdx++]);
          string field(arg[argIdx++]);
          pair <string,FieldName> pair_name(name,string_to_field(field));
          if (++argIdx < narg) { // keyword average
            esetData_[pair_name] = ELEMENTSET_AVERAGE;
          }
          else {
            esetData_[pair_name] = ELEMENTSET_TOTAL;
          }
          match = true;
        }

        else if (strcmp(arg[argIdx],"now")==0) {
          argIdx++;
          double dt = 1.0;
          if (argIdx < narg) {
            dt = atof(arg[argIdx++]);
          }
          update_time(dt);
          update_step();
          outputNow_ = true;
          this->output();
          outputNow_ = false;
          match = true;
        }
        else
          if (strcmp(arg[argIdx],"index")==0) {
            argIdx++;
            if (strcmp(arg[argIdx],"step")==0) { outputTime_ = false; }
            else                               { outputTime_ = true; }
          match = true;
        }
        else {
          outputPrefix_ = arg[argIdx++];
          outputFrequency_ = atoi(arg[argIdx++]);
          bool ensight_output = false, full_text_output = false;
          bool text_output = false, vect_comp = false, tensor_comp = false;
          int rank = lammpsInterface_->comm_rank();
          for (int i = argIdx; i<narg; ++i) {
            if      (strcmp(arg[i],"full_text")==0) full_text_output = true;
            else if (strcmp(arg[i],"text")==0)           text_output = true;
            else if (strcmp(arg[i],"binary")==0)      ensight_output = true;
            else if (strcmp(arg[i],"vector_components")==0) vect_comp = true;
            else if (strcmp(arg[i],"tensor_components")==0) tensor_comp = true;
            else { throw ATC_Error(" output: unknown keyword ");  }
          }
          if (outputFrequency_>0) {
            set<int> otypes;
            if (full_text_output || text_output) {
              lammpsInterface_->print_msg_once("Warning : text output can create _LARGE_ files");
            }
            if (full_text_output) otypes.insert(FULL_GNUPLOT);
            if (text_output)      otypes.insert(GNUPLOT);
            if (ensight_output)   otypes.insert(ENSIGHT);
            if (ntracked() > 0) {
               string fstem = field_to_string(SPECIES_CONCENTRATION);
               string istem = field_to_intrinsic_name(SPECIES_CONCENTRATION);
               vector<string> tnames = tracked_names();
               vector<string> fnames;
               vector<string> inames;
               for (unsigned int i = 0; i < tnames.size(); i++) {
                 fnames.push_back(fstem+tnames[i]);
                 inames.push_back(istem+tnames[i]);
               }
               feEngine_->add_field_names(fstem,fnames);
               feEngine_->add_field_names(istem,inames);
            }
            feEngine_->initialize_output(rank,outputPrefix_,otypes);
            if (vect_comp)
              feEngine_->output_manager()
                ->set_option(OUTPUT_VECTOR_COMPONENTS,true);
            if (tensor_comp)
              feEngine_->output_manager()
                ->set_option(OUTPUT_TENSOR_COMPONENTS,true);
          }
          match = true;
        }
      }
    else if (strcmp(arg[argIdx],"write")==0) {
      argIdx++;
      FieldName thisField;
      int thisIndex;
      parse_field(arg,argIdx,thisField,thisIndex);
      string nsetName(arg[argIdx++]);
      string filename(arg[argIdx++]);
      stringstream  f;
      set<int> nodeSet = (feEngine_->fe_mesh())->nodeset(nsetName);
      set<int>::const_iterator iset;
      const DENS_MAT & field =(fields_.find(thisField)->second).quantity();
      for (iset = nodeSet.begin(); iset != nodeSet.end(); iset++) {
        int inode = *iset;
        double v = field(inode,thisIndex);
        f << inode << " " << std::setprecision(17) << v << "\n";
      }
      LammpsInterface::instance()->write_file(filename,f.str());
      match = true;
    }
    // add a species for tracking
    /*! \page man_add_species fix_modify AtC add_species
      \section syntax
      fix_modify_AtC add_species <TAG> <group|type> <ID> \n
      - <TAG> = tag for tracking a species \n
      - group|type = LAMMPS defined group or type of atoms \n
      - <ID> = name of group or type number \n
      \section examples
      <TT> fix_modify AtC add_species gold type 1 </TT> \n
      <TT> group GOLDGROUP type 1 </TT> \n
      <TT> fix_modify AtC add_species gold group GOLDGROUP </TT>
      \section description
      Associates a tag with all atoms of a specified type or within a specified group. \n
      \section restrictions
      \section related
      \section default
      No defaults for this command.
    */
    else if (strcmp(arg[argIdx],"add_species")==0) {
      argIdx++;
      string speciesTag = arg[argIdx];
      string tag = arg[argIdx];
      argIdx++;
      if (strcmp(arg[argIdx],"group")==0) {
        if (narg-argIdx == 2) {
          string name = arg[++argIdx];
          int id = lammpsInterface_->group_bit(name);
          groupList_.push_back(id);
          groupNames_.push_back(tag);
        }
        else {
          while (++argIdx < narg) {
            string name = arg[argIdx];
            int id = lammpsInterface_->group_bit(name);
            string tag = speciesTag+"-"+name;
            groupList_.push_back(id);
            groupNames_.push_back(tag);
          }
        }
      }
      else if (strcmp(arg[argIdx],"type")==0) {
        if (narg-argIdx == 2) {
          int id = atoi(arg[++argIdx]);
          typeList_.push_back(id);
          typeNames_.push_back(tag);
        }
        else {
          while (++argIdx < narg) {
            int id = atoi(arg[argIdx]);
            string tag = speciesTag+"_"+to_string(id);
            typeList_.push_back(id);
            typeNames_.push_back(tag);
          }
        }
      }
      else {
        throw ATC_Error("ATC_Method: add_species only handles groups or types"); }
      match = true;
    }

    // remove species from tracking

    /*! \page man_remove_species fix_modify AtC remove_species
      \section syntax
      fix_modify_AtC delete_species <TAG> \n

      - <TAG> = tag for tracking a species \n
      \section examples
      <TT> fix_modify AtC remove_species gold </TT> \n
      \section description
      Removes tag designated for tracking a specified species. \n
      \section restrictions
      \section related
      \section default
      No defaults for this command.
    */
    else if (strcmp(arg[argIdx],"delete_species")==0) {
      argIdx++;
      string tag = arg[argIdx++];
      if (strcmp(arg[argIdx],"group")==0) {
        for (unsigned int j = 0; j < groupList_.size(); j++) {
          if (tag == groupNames_[j]) {
            groupList_.erase(groupList_.begin()+j);
            groupNames_.erase(groupNames_.begin()+j);
            break;
          }
        }
      }
      else if (strcmp(arg[argIdx],"type")==0) {
        for (unsigned int j = 0; j < typeList_.size(); j++) {
          if (tag == typeNames_[j]) {
            typeList_.erase(typeList_.begin()+j);
            typeNames_.erase(typeNames_.begin()+j);
            break;
          }
        }
      }
      else {
        throw ATC_Error("ATC_Method: delete_species only handles groups or types"); }
      match = true;

    }

    // add a molecule for tracking
    /*! \page man_add_molecule fix_modify AtC add_molecule
      \section syntax
      fix_modify_AtC add_molecule <small|large> <TAG> <GROUP_NAME> \n

      - small|large = can be small if molecule size < cutoff radius, must be large otherwise \n
      - <TAG> = tag for tracking a species \n
      - <GROUP_NAME> = name of group that tracking will be applied to \n
      \section examples
      <TT> group WATERGROUP type 1 2 </TT> \n
      <TT> fix_modify AtC add_molecule small water WATERGROUP </TT> \n
      \section description
      Associates a tag with all molecules corresponding to a specified group. \n
      \section restrictions
      \section related
      \section default
      No defaults for this command.
    */
    else if (strcmp(arg[argIdx],"add_molecule")==0) {
      argIdx++;
      MolSize size;
      if (strcmp(arg[argIdx],"small")==0) {
        size = MOL_SMALL;
        //needXrefProcessorGhosts_ = true;
        needProcGhost_ = true;
      }
      else
        throw ATC_Error("ATC_CouplingMass:  Bad molecule size in add_molecule");
      argIdx++;
      string moleculeTag = arg[argIdx];

      argIdx++;
      int groupBit = lammpsInterface_->group_bit(arg[argIdx]);
      moleculeIds_[moleculeTag] = pair<MolSize,int>(size,groupBit);
      match = true;
    }

    // remove molecule from tracking
    /*! \page man_remove_molecule fix_modify AtC remove_molecule
      \section syntax
      fix_modify_AtC remove_molecule <TAG> \n

      - <TAG> = tag for tracking a molecule type \n
      \section examples
      <TT> fix_modify AtC remove_molecule water </TT> \n
      \section description
      Removes tag designated for tracking a specified set of molecules. \n
      \section restrictions
      \section related
      \section default
      No defaults for this command.
    */
    else if (strcmp(arg[argIdx],"remove_molecule")==0) {
      argIdx++;
      string moleculeTag = arg[argIdx];
      moleculeIds_.erase(moleculeTag);

      taggedDensMan_.erase(moleculeTag);
    }

      /*! \page man_boundary fix_modify AtC boundary
        \section syntax
        fix_modify AtC boundary type <atom-type-id>
        - <atom-type-id> = type id for atoms that represent a fictitious
        boundary internal to the FE mesh
        \section examples
        <TT> fix_modify AtC boundary type ghost_atoms </TT>
        \section description
        Command to define the atoms that represent the fictitious
        boundary internal to the FE mesh. For fully overlapped MD/FE
        domains with periodic boundary conditions no boundary atoms should
        be defined.
        \section restrictions
        \section default
        none
      */
      else if (strcmp(arg[argIdx],"boundary")==0) {
        argIdx++;
        groupTagGhost_ = arg[argIdx++];
        match = true;
      }

      /*! \page man_internal_atom_integrate fix_modify AtC internal_atom_integrate
        \section syntax
        fix_modify AtC internal_atom_integrate <on | off>
        <TT> fix_modify AtC internal_atom_integrate on </TT>
        \section description
        Has AtC perform time integration for the atoms in the group on which it operates.  This does not include boundary atoms.
        \section restrictions
        AtC must be created before any fixes doing time integration.  It must be on for coupling methods which impose constraints on velocities during the first verlet step, e.g. control momentum glc_velocity.
        \section default
        on for coupling methods, off for post-processors
        off
       */
      else if (strcmp(arg[argIdx],"internal_atom_integrate")==0) {
        argIdx++;
        if (strcmp(arg[argIdx],"off")==0) {
          integrateInternalAtoms_ = false;
          match = true;
        }
        else {
          integrateInternalAtoms_ = true;
          match = true;
        }
      }

      /*! \page man_internal_element_set fix_modify AtC internal_element_set
        \section syntax
        fix_modify AtC internal_element_set <element-set-name>
        - <element-set-name> = name of element set defining internal region, or off
        \section examples
        <TT> fix_modify AtC internal_element_set myElementSet </TT>
        <TT> fix_modify AtC internal_element_set off </TT>
        \section description
        Enables AtC to base the region for internal atoms to be an element set.
        If no ghost atoms are used, all the AtC atoms must be constrained to remain
        in this element set by the user, e.g., with walls.  If boundary atoms are
        used in conjunction with Eulerian atom maps
        AtC will partition all atoms of a boundary or internal type to be of type internal
        if they are in the internal region or to be of type boundary otherwise.
        \section restrictions
        If boundary atoms are used in conjunction with Eulerian atom maps, the Eulerian
        reset frequency must be an integer multiple of the Lammps reneighbor frequency
        \section related
        see \ref atom_element_map_type and \ref boundary
        \section default
        off
       */
      else if (strcmp(arg[argIdx],"internal_element_set")==0) {
        argIdx++;
        if (strcmp(arg[argIdx],"off")==0) {
          internalElementSet_ = "";
          match = true;
        }
        else {
          internalElementSet_ = string(arg[argIdx]);
          const set<int> & elementSet((feEngine_->fe_mesh())->elementset(internalElementSet_)); // check it exists and is not trivial
          if (elementSet.size()==0) throw ATC_Error("internal_element_set - element set " + internalElementSet_ + " has no elements");
          match = true;
        }
      }

    /*! \page man_atom_weight fix_modify AtC atom_weight
      \section syntax
      fix_modify AtC atom_weight <method> <arguments>
        - <method> = \n
          value: atoms in specified group assigned constant value given \n
          lattice: volume per atom for specified lattice type (e.g. fcc) and parameter \n
          element: element volume divided among atoms within element \n
          region: volume per atom determined based on the atom count in the MD regions and their volumes. Note: meaningful only if atoms completely fill all the regions. \n
          group: volume per atom determined based on the atom count in a group and its volume\n
          read_in: list of values for atoms are read-in from specified file \n
      \section examples
       <TT> fix_modify atc atom_weight constant myatoms 11.8 </TT> \n
       <TT> fix_modify atc atom_weight lattice </TT> \n
       <TT> fix_modify atc atom_weight read-in atm_wt_file.txt </TT> \n
      \section description
       Command for assigning the value of atomic weights used for atomic integration in
       atom-continuum coupled simulations.
      \section restrictions
      Use of lattice option requires a lattice type and parameter is already specified.
      \section related
      \section default
      lattice
    */
      else if (strcmp(arg[argIdx],"atom_weight")==0) {
        argIdx++;
        if (strcmp(arg[argIdx],"constant")==0) {
          argIdx++;
          atomWeightType_ = USER;
          int groupbit = -1;
          if (strcmp(arg[argIdx],"all")==0) {
          }
          else {
            groupbit = lammpsInterface_->group_bit(arg[argIdx]);
          }
          argIdx++;
          double value = atof(arg[argIdx]);
          Valpha_[groupbit] = value;
          match = true;
        }
        else if (strcmp(arg[argIdx],"lattice")==0) {
          atomWeightType_ = LATTICE;
          match = true;
        }
        else if (strcmp(arg[argIdx],"element")==0) {
          atomWeightType_ = ELEMENT;
          match = true;
        }
        else if (strcmp(arg[argIdx],"region")==0) {
          atomWeightType_ = REGION;
          match = true;
        }
        else if (strcmp(arg[argIdx],"group")==0) {
          atomWeightType_ = GROUP;
          match = true;
        }
        else if (strcmp(arg[argIdx],"multiscale")==0) {
          atomWeightType_ = MULTISCALE;
          match = true;
        }
        else if (strcmp(arg[argIdx],"node")==0) {
          atomWeightType_ = NODE;
          match = true;
        }
        else if (strcmp(arg[argIdx],"node_element")==0) {
          atomWeightType_ = NODE_ELEMENT;
          match = true;
        }
        else if (strcmp(arg[argIdx],"read_in")==0) {
          atomWeightType_ = READ_IN;
          argIdx++;
          atomicWeightsFile_ = arg[argIdx];
          match = true;
        }
        if (match) {
          needReset_ = true;
        }
      }

    /*! \page man_decomposition fix_modify AtC decomposition
      \section syntax
      fix_modify AtC decomposition <type>
        - <type> = \n
          replicated_memory: nodal information replicated on each processor \n
          distributed_memory: only owned nodal information on processor  \n
      \section examples
       <TT> fix_modify atc decomposition distributed_memory </TT> \n
      \section description
       Command for assigning the distribution of work and memory for parallel runs.
      \section restrictions
      replicated_memory is appropriate for simulations were the number of nodes << number of atoms
      \section related
      \section default
      replicated_memory
    */
      else if (strcmp(arg[argIdx],"decomposition")==0) {
        argIdx++;
        if (strcmp(arg[argIdx],"replicated_memory")==0) {
          domainDecomposition_ = REPLICATED_MEMORY;
          match = true;
        }
        else if (strcmp(arg[argIdx],"distributed_memory")==0) {
          domainDecomposition_ = DISTRIBUTED_MEMORY;
          match = true;
        }
      }

    /*! \page man_write_atom_weights fix_modify AtC write_atom_weights
      \section syntax
      fix_modify AtC write_atom_weights <filename> <frequency>
        - <filename> = name of file that atomic weights are written to \n
        - <frequency> = how often writes will occur \n
      \section examples
       <TT> fix_modify atc write_atom_weights atm_wt_file.txt 10 </TT> \n
      \section description
       Command for writing the values of atomic weights to a specified file.
      \section restrictions
      \section related
      \section default
    */
      else if (strcmp(arg[argIdx],"write_atom_weights")==0) {
        argIdx++;
        atomicWeightsFile_ = arg[argIdx];
        argIdx++;
        atomicWeightsWriteFrequency_ = atoi(arg[argIdx]);
        atomicWeightsWriteFlag_ = true;
        match = true;
      }


      /*! \page man_reset_time fix_modify AtC reset_time
      \section syntax
      fix_modify AtC reset_time <value>
      \section examples
       <TT> fix_modify atc reset_time 0.0 </TT> \n
      \section description
      Resets the simulation time counter.
      \section restrictions
      \section related
      \section default
      */
      else if (strcmp(arg[argIdx],"reset_time")==0) {
        argIdx++;
        set_time();
        if (narg > argIdx) {
          double time = atof(arg[argIdx]);
          set_time(time);
        }
        match = true;
      }

      /*! \page man_reset_time fix_modify AtC kernel_bandwidth
      \section syntax
      fix_modify AtC kernel_bandwidth <value>
      \section examples
       <TT> fix_modify atc kernel_bandwidth 8 </TT> \n
      \section description
      Sets a maximum parallel bandwidth for the kernel functions during parallel communication.  If the command is not issued, the default will be to assume the bandwidth of the kernel matrix corresponds to the number of sampling locations.
      \section restrictions
      Only is used if kernel functions are being used.
      \section related
      \section default
      Number of sample locations.
      */
      else if (strcmp(arg[argIdx],"kernel_bandwidth")==0) {
        argIdx++;
        accumulantBandwidth_ = atoi(arg[argIdx]);
        match = true;
      }

      /*! \page man_reset_atomic_reference_positions fix_modify AtC reset_atomic_reference_positions
      \section syntax
      fix_modify AtC reset_atomic_reference_positions
      \section examples
       <TT> fix_modify atc reset_atomic_reference_positions
      \section description
      Resets the atomic positions ATC uses to perform point to field operations.
      In can be used to use perfect lattice sites in ATC but a thermalized or
      deformed lattice in LAMMPS.
      \section restrictions
      \section related
      \section default
      Default is off
      */
      else if (strcmp(arg[argIdx],"reset_atomic_reference_positions")==0) {
        argIdx++;
        xRefFile_ = arg[argIdx];
        readXref_ = true;
        match = true;
      }

      /*! \page man_set fix_modify AtC set
        \section syntax
        fix_modify AtC set reference_potential_energy <value_or_filename(optional)>
        - value (double) : optional user specified zero point for PE in native LAMMPS energy units \n
        - filename (string) : optional user specified string for file of nodal PE values to be read-in
        \section examples
        <TT> fix_modify AtC set reference_potential_energy </TT> \n
        <TT> fix_modify AtC set reference_potential_energy -0.05 </TT> \n
        <TT> fix_modify AtC set reference_potential_energy myPEvalues </TT> \n
        \section description
        Used to set various quantities for the post-processing algorithms.
        It sets the zero point for the potential energy density using
        the value provided for all nodes, or from the current
        configuration of the lattice if no value is provided, or
        values provided within the specified filename.
        \section restrictions
        Must be used with the hardy/field type of AtC fix
        ( see \ref man_fix_atc )
        \section related
        \section default
        Defaults to lammps zero point i.e. isolated atoms
      */
      else if (strcmp(arg[argIdx],"set")==0) {
        argIdx++;
        if (strcmp(arg[argIdx],"reference_potential_energy")==0) {
          argIdx++;
          setRefPE_ = true;
          if (narg > argIdx) {
            string a(arg[argIdx]);
            if (is_dbl(a)) {
              double value = atof(arg[argIdx]);
              refPEvalue_ = value;
              setRefPEvalue_ = true;
            }
            else {
              nodalRefPEfile_ = arg[argIdx];
              readRefPE_ = true;
            }
          }
          match = true;
        }
      } // end "set"



      /*! \page man_atom_element_map fix_modify AtC atom_element_map
      \section syntax
      fix_modify AtC atom_element_map  <eulerian|lagrangian> <frequency> \n
      - frequency (int) : frequency of updating atom-to-continuum maps based on the
      current configuration - only for eulerian
      \section examples
      <TT> fix_modify atc atom_element_map eulerian 100 </TT>
      \section description
      Changes frame of reference from eulerian to lagrangian and sets the
      frequency for which the map from atoms to elements is reformed and
      all the attendant data is recalculated.
      \section restrictions
      Cannot change map type after initialization.
      \section related
      \section default
      lagrangian
      */
      else if (strcmp(arg[argIdx],"atom_element_map")==0) {
        argIdx++;
        if (strcmp(arg[argIdx],"eulerian")==0) {
          atomToElementMapType_ = EULERIAN;
          argIdx++;
          atomToElementMapFrequency_ = atoi(arg[argIdx]);
        }
        else {
          atomToElementMapType_ = LAGRANGIAN;
          atomToElementMapFrequency_ = 0;
        }
        match = true;
        needReset_ = true;
      }

      /*! \page man_read_restart fix_modify AtC read_restart
      \section syntax
      fix_modify AtC read_restart [file_name]  \n
      \section examples
      <TT> fix_modify AtC read_restart ATC_state </TT> \n
      \section description
      Reads the current state of the fields from a named text-based restart
      file.
      \section restrictions
      The restart file only contains fields and their time derivatives.
      The reference positions of the atoms and the commands that initialize
      the fix are not saved e.g. an identical mesh containing the same atoms
      will have to be recreated.
      \section related
      see write_restart \ref man_write_restart
      \section default
      none
      */
      else if (strcmp(arg[argIdx],"read_restart")==0) {
        argIdx++;
        restartFileName_  = arg[argIdx];
        useRestart_ = true;
        match = true;
      }

      /*! \page man_write_restart fix_modify AtC write_restart
      \section syntax
      fix_modify AtC write_restart [file_name]  \n
      \section examples
      <TT> fix_modify AtC write_restart restart.mydata </TT> \n
      \section description
      Dumps the current state of the fields to a named text-based restart file.
      This done when the command is invoked and not repeated, unlike the
      similar lammps command.
      \section restrictions
      The restart file only contains fields and their time derivatives.
      The reference positions of the atoms and the commands that initialize
      the fix are not saved e.g. an identical mesh containing the same atoms
      will have to be recreated.
      \section related
      see read_restart \ref man_read_restart
      \section default
      none
      */
      else if (strcmp(arg[argIdx],"write_restart")==0) {
        argIdx++;
        string restartFileName(arg[argIdx]);
        RESTART_LIST data;
        write_restart_data(restartFileName,data);
        match = true;
      }

    } // end else

    return match; // return to FixATC

  }

  //--------------------------------------------------
  // helper function for parser
  // handles : "displacement x"  and "temperature" by indexing argIdx
  // for fluxes : only normal fluxes can be prescribed
  //--------------------------------------------------
  void ATC_Method::parse_field(/*const*/ char ** args, int & argIdx,
                               FieldName & thisField)
  {
    string thisName = args[argIdx++];
    thisField = string_to_field(thisName);
    map<FieldName,int>::const_iterator iter = fieldSizes_.find(thisField);
    if (iter == fieldSizes_.end()) {
      throw ATC_Error("Bad field name: "+thisName);
    }
  }

  //--------------------------------------------------
  // helper function for parser
  // handles : "displacement x"  and "temperature" by indexing argIdx
  // for fluxes : only normal fluxes can be prescribed
  //--------------------------------------------------
  void ATC_Method::parse_field(/*const*/ char ** args, int & argIdx,
                           FieldName & thisField, int & thisIndex)
  {
    string thisName = args[argIdx++];
    if (args[argIdx] == nullptr) {
      throw ATC_Error("Need to give field '"+thisName+"' more args");
    }
    thisField = string_to_field(thisName);
    map<FieldName,int>::const_iterator iter = fieldSizes_.find(thisField);
    if (iter == fieldSizes_.end()) {
      throw ATC_Error("Bad field name: "+thisName);
    }
    string thisDim = args[argIdx];
    thisIndex = 0;
    if (string_to_index(thisDim,thisIndex)) {
      if ( !( thisIndex < fieldSizes_[thisField]) ) {
        throw ATC_Error("Bad field dimension "+thisDim);
      }
      argIdx++;
    }
  }

  //-------------------------------------------------------------------
  // this sets the peratom output

  void ATC_Method::update_peratom_output()
  {
    perAtomArray_ = perAtomOutput_;
    // copy values
    for (int i = 0; i < lammpsInterface_->nlocal(); i++) {
      for (int j = 0; j < nsd_; j++) {
        perAtomOutput_[i][j] = xref_[i][j];
      }
      for (int j = nsd_; j < sizePerAtomCols_; j++) {
        perAtomOutput_[i][j] = 0;
      }
    }
    int indx = nsd_;
    if (atomVolume_->nRows() > 0) { // kernel Hardy does not compute these
      const DIAG_MAT & myAtomicWeights(atomVolume_->quantity());
      for (int i = 0; i < nLocal_; i++) {
        double wg = myAtomicWeights(i,i);
        if (wg > 0) {
          int ii = internalToAtom_(i);
          perAtomOutput_[ii][indx] = 1./wg;
        }
      }
    }
  }

  void ATC_Method::adjust_xref_pbc()
  {

    int nlocal = lammpsInterface_->nlocal();
    int xperiodic = lammpsInterface_->xperiodic();
    int yperiodic = lammpsInterface_->yperiodic();
    int zperiodic = lammpsInterface_->zperiodic();
    double **x = lammpsInterface_->xatom();
    double boxxlo,boxxhi;
    double boxylo,boxyhi;
    double boxzlo,boxzhi;

    lammpsInterface_->box_bounds(boxxlo,boxxhi,
                                 boxylo,boxyhi,
                                 boxzlo,boxzhi);
//  bool changed = false;
    for (int i = 0; i < nlocal; i++) {
      if (xperiodic) {
        if (x[i][0] < boxxlo) {
          xref_[i][0] += Xprd_;
//        changed = true;
        }
        if (x[i][0] >= boxxhi) {
          xref_[i][0] -= Xprd_;
//        changed = true;
        }
      }

      if (yperiodic) {
        if (x[i][1] < boxylo) {
          xref_[i][1] += Yprd_;
//        changed = true;
        }
        if (x[i][1] >= boxyhi) {
          xref_[i][1] -= Yprd_;
//        changed = true;
        }
      }

      if (zperiodic) {
        if (x[i][2] < boxzlo) {
          xref_[i][2] += Zprd_;
//        changed = true;
        }
        if (x[i][2] >= boxzhi) {
          xref_[i][2] -= Zprd_;
//        changed = true;
        }
      }
    }

    // propagate reset if needed
    if (atomToElementMapType_ == LAGRANGIAN) {
      if (atomCoarseGrainingPositions_) {
        atomCoarseGrainingPositions_->force_reset();
      }
    }
    else if (atomReferencePositions_) {
      atomReferencePositions_->force_reset();
    }

  }
  //-------------------------------------------------------------------
  void ATC_Method::initialize()
  {
    feEngine_->partition_mesh();
    // initialized_ is set to true by derived class initialize()
    // localStep_ is a counter within a run or minimize
    localStep_ = 0;
    // STEP 1)  get basic information data from Lammps/fix
    // 1a)  group ids for all internal atoms
    groupbit_ = lammpsInterface_->group_bit(groupTag_);

    // 1b) group ids for ghost atoms
    groupbitGhost_ = 0;
    if (!groupTagGhost_.empty()) {
      groupbitGhost_ = lammpsInterface_->group_bit(groupTagGhost_);
    }

    // 1c) periodicity and box bounds/lengths
    if (!initialized_) {

      lammpsInterface_->box_periodicity(periodicity[0],
                                            periodicity[1],
                                            periodicity[2]);
      lammpsInterface_->box_bounds(box_bounds[0][0],box_bounds[1][0],
                                       box_bounds[0][1],box_bounds[1][1],
                                       box_bounds[0][2],box_bounds[1][2]);
      for (int k = 0; k < nsd_; k++) {
        box_length[k] = box_bounds[1][k] - box_bounds[0][k];
      }

      lammpsInterface_->set_reference_box();

      // get periodicity data from lammps for parallel exchange to adjust for periodicity
      Xprd_ = lammpsInterface_->domain_xprd();
      Yprd_ = lammpsInterface_->domain_yprd();
      Zprd_ = lammpsInterface_->domain_zprd();
//    box_length[0] = Xprd_;
//    box_length[1] = Yprd_;
//    box_length[2] = Zprd_;
      XY_ = lammpsInterface_->domain_xy();
      XZ_ = lammpsInterface_->domain_xz();
      YZ_ = lammpsInterface_->domain_yz();
    }

    // STEP 2 computational geometry
    // 2a) get basic information from continuum/FE
    this->set_continuum_data();

    // STEP 2b) set up data structures for computational geometry
    if (this->reset_methods()) {
      // clear memory manager
      interscaleManager_.clear_temporary_data();
      atomVolume_ = nullptr;

      // reference positions and energy
      if (!initialized_) {
        double **x = lammpsInterface_->xatom();
        for (int i = 0; i < lammpsInterface_->nmax(); i++) {
          for (int j = 0; j < nsd_; j++) {
            xref_[i][j] = x[i][j];
          }
        }

        // re-write non-ghosts' xref with values from a file
        if (readXref_) {
          bool success = read_atomic_ref_positions(xRefFile_.c_str());
          if (!success)
            throw ATC_Error("Error reading atomic reference positions");
          readXref_ = false;
        }

        // ensure initial configuration is consistent with element set

        if (internalElementSet_.size() && groupbitGhost_) {
          int *mask = lammpsInterface_->atom_mask();
          int nlocal = lammpsInterface_->nlocal();
          const FE_Mesh * feMesh = feEngine_->fe_mesh();
          const set<int> & elementSet(feMesh->elementset(internalElementSet_));
          int element;
          DENS_VEC coords(nsd_);
          for (int i = 0; i < nlocal; ++i) {
            if (mask[i] & groupbit_ || mask[i] & groupbitGhost_) {
              for (int j = 0; j < nsd_; j++) {
                coords(j) = xref_[i][j];
              }
              element = feMesh->map_to_element(coords);
              if (elementSet.find(element) == elementSet.end()) {
                mask[i] |= groupbitGhost_;
                mask[i] &= ~groupbit_;
              }
              else {
                mask[i] &= ~groupbitGhost_;
                mask[i] |= groupbit_;
              }
            }
          }
        }

        // set up maps from lammps to atc indexing
        reset_nlocal();
      }

      this->set_computational_geometry();
    }

    // 2c) basic data regarding atomic system, e.g. atom coordinates
    if (atomToElementMapType_ == EULERIAN) {
      reset_coordinates();
    }

    // STEP 3) set up variables which will be integrated in time
    this->construct_time_integration_data();

    // STEP 4) instantiate all the various specific algorithms and methods
    this->construct_methods();

    // STEP 5) construct dependency-managed data
    // 5b) all other transfer operators
    // needs to be done before every run in case options have changed or the atoms have been changed by the user

    if (this->reset_methods()) {
      // construct all the needed data structures
      this->construct_transfers();

      // allocate all space needed for lammps arrays
      interscaleManager_.grow_arrays(lammpsInterface_->nmax());
    }
    // reset all computes invoked flags and lammps data
    interscaleManager_.lammps_force_reset();

    // STEP 6) initialize data
    // 6b) size quantities which use pack_comm
    interscaleManager_.size_comm_quantities();

    // 6c) set coarse-graining functions and atomic weights
    if (!initialized_) {
      // FE_Engine allocates all required memory
      // assume initial atomic position is the reference position for now

      // \int_\Omega N_I dV : static if the mesh is
      NodeVolumes_.reset(nNodes_,nNodes_);
      invNodeVolumes_.reset(nNodes_,nNodes_);
      feEngine_->compute_lumped_mass_matrix(NodeVolumes_);
      invNodeVolumes_.set_quantity() = NodeVolumes_.inv();
    }
    atomVolume_->set_reset();

    // 6d) reference values
    this->set_reference_potential_energy();

    // 6e) atomic output for 0th step
    update_peratom_output();

    massMatInv_.reset(nNodes_,nNodes_);
    feEngine_->compute_lumped_mass_matrix(massMatInv_);
    for (int i = 0; i < nNodes_; ++i)  {
      massMatInv_(i,i) = 1./massMatInv_(i,i);
    }

    // clear need for resets
    needReset_ = false;

  }
  //-------------------------------------------------------------------
  void ATC_Method::set_continuum_data()
  {
    // initialize finite element engine and get basic properties
    if (!initialized_) {
      feEngine_->initialize();
      if (nsd_!=feEngine_->nsd()) {
        throw ATC_Error("Spatial dimensions inconsistent between LAMMPS and ATC");
      }
      nNodes_ = feEngine_->num_nodes();
    }
  }

  //-------------------------------------------------------------------
  void ATC_Method::set_computational_geometry()
  {
    // set positions used for coarse-graining operators




    if (!initialized_) {
      if (atomToElementMapType_ == EULERIAN) {
        FundamentalAtomQuantity * atomPositionsAll = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION,ALL);
        ClonedAtomQuantity<double> * myAtomPositions =
          new ClonedAtomQuantity<double>(this,atomPositionsAll,INTERNAL);
        atomCoarseGrainingPositions_ = myAtomPositions;
        interscaleManager_.add_per_atom_quantity(myAtomPositions,
                                                 "AtomicCoarseGrainingPositions");

        if (trackDisplacement_) {
          XrefWrapper * myAtomReferencePositions = new XrefWrapper(this);
          atomReferencePositions_ = myAtomReferencePositions;
          interscaleManager_.add_per_atom_quantity(myAtomReferencePositions,
                                                   "AtomicReferencePositions");
          atomReferencePositions_->set_memory_type(PERSISTENT);
        }

        if (groupbitGhost_) {
          myAtomPositions = new ClonedAtomQuantity<double>(this,atomPositionsAll,GHOST);
          atomGhostCoarseGrainingPositions_ = myAtomPositions;
          interscaleManager_.add_per_atom_quantity(myAtomPositions,
                                                   "AtomicGhostCoarseGrainingPositions");
        }
        if(needProcGhost_){
          FundamentalAtomQuantity * atomPositionsAll = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION,PROC_GHOST);
          ClonedAtomQuantity<double> * myAtomPositions =
            new ClonedAtomQuantity<double>(this,atomPositionsAll,PROC_GHOST);
          atomProcGhostCoarseGrainingPositions_ = myAtomPositions;
          interscaleManager_.add_per_atom_quantity(myAtomPositions,
                                                   "AtomicProcGhostCoarseGrainingPositions");
        }
      }
      else {
        XrefWrapper * myAtomPositions = new XrefWrapper(this);
        atomCoarseGrainingPositions_ = myAtomPositions;
        interscaleManager_.add_per_atom_quantity(myAtomPositions,
                                                 "AtomicCoarseGrainingPositions");
        atomReferencePositions_ = atomCoarseGrainingPositions_;

        if (groupbitGhost_) {
          myAtomPositions = new XrefWrapper(this,GHOST);
          atomGhostCoarseGrainingPositions_ = myAtomPositions;
          interscaleManager_.add_per_atom_quantity(myAtomPositions,
                                                   "AtomicGhostCoarseGrainingPositions");
        }
        if (needProcGhost_) {
          XrefWrapper * myAtomPositions = new XrefWrapper(this);
          atomProcGhostCoarseGrainingPositions_ = myAtomPositions;
          interscaleManager_.add_per_atom_quantity(myAtomPositions,
                                                   "AtomicProcGhostCoarseGrainingPositions");
        }
      }
      atomCoarseGrainingPositions_->set_memory_type(PERSISTENT);
      if (atomGhostCoarseGrainingPositions_) atomGhostCoarseGrainingPositions_->set_memory_type(PERSISTENT);
      if (atomProcGhostCoarseGrainingPositions_) atomProcGhostCoarseGrainingPositions_->set_memory_type(PERSISTENT);
    }

    // Add in atom to element map if using shape functions
    if (needsAtomToElementMap_) {
      atomElement_ = new AtomToElementMap(this);
      interscaleManager_.add_per_atom_int_quantity(atomElement_,"AtomElement");
    }
  }

  //-------------------------------------------------------------------
  void ATC_Method::construct_methods()
  {

    if (this->reset_methods()) {
      if (atomTimeIntegrator_) delete atomTimeIntegrator_;
      if (integrateInternalAtoms_) {
        atomTimeIntegrator_ = new AtomTimeIntegratorType(this,INTERNAL);
      }
      else {
        atomTimeIntegrator_ = new AtomTimeIntegrator();
      }

      // set up integration schemes for ghosts
      ghostManager_.construct_methods();
    }
  }

  //-------------------------------------------------------------------
  void ATC_Method::construct_transfers()
  {
    this->construct_interpolant();

    this->construct_molecule_transfers();

    atomTimeIntegrator_->construct_transfers();
    ghostManager_.construct_transfers();



  }
  //-------------------------------------------------------------------
  PerAtomDiagonalMatrix<double> * ATC_Method::create_atom_volume()
  {
    if (atomVolume_) {
      return atomVolume_;
    }
    else {
      // set variables to compute atomic weights
      DENS_MAN * nodalVolume(nullptr);
      switch (atomWeightType_) {
      case USER:
        atomVolume_ = new AtomVolumeUser(this,Valpha_);
        break;
      case LATTICE:
        atomVolume_ = new AtomVolumeLattice(this);
        break;
      case ELEMENT:
        atomVolume_ = new AtomVolumeElement(this);
        break;
      case REGION:
        atomVolume_ = new AtomVolumeRegion(this);
        break;
      case GROUP:
        atomVolume_ = new AtomVolumeGroup(this,Valpha_);
        break;
      case MULTISCALE:
        if (!shpFcn_) {
          throw ATC_Error("ATC_Method::create_atom_volume - Multiscale algorithm requires an interpolant");
        }
        nodalVolume = new NodalAtomVolume(this,shpFcn_);
        interscaleManager_.add_dense_matrix(nodalVolume,"NodalAtomVolume");
        atomVolume_ = new FtaShapeFunctionProlongationDiagonalMatrix(this,nodalVolume,shpFcn_);
        break;
      case NODE:
        if (!shpFcn_) {
          throw ATC_Error("ATC_Method::create_atom_volume - Node algorithm requires an interpolant");
        }
        nodalVolume = new NodalVolume(this,shpFcn_);
        interscaleManager_.add_dense_matrix(nodalVolume,"NodalVolume");
        atomVolume_ = new FtaShapeFunctionProlongationDiagonalMatrix(this,nodalVolume,shpFcn_);
        break;
      case NODE_ELEMENT:
        if (!shpFcn_) {
          throw ATC_Error("ATC_Method::create_atom_volume - Node-Element algorithm requires an interpolant");
        }
        nodalVolume = new NodalAtomVolumeElement(this,shpFcn_);
        interscaleManager_.add_dense_matrix(nodalVolume,"NodalAtomVolumeElement");
        atomVolume_ = new FtaShapeFunctionProlongationDiagonalMatrix(this,nodalVolume,shpFcn_);
        break;
      case READ_IN:
        atomVolume_ = new AtomVolumeFile(this,atomicWeightsFile_);
        break;
      }
      if (atomVolume_) {
        interscaleManager_.add_per_atom_diagonal_matrix(atomVolume_,"AtomVolume");
      }
      else {
        throw ATC_Error("ATC_Method::create_atom_volume - bad option for atom volume algorithm");
      }

      return atomVolume_;
    }
  }
  //--------------------------------------------------------
  void ATC_Method::init_integrate()
  {
    atomTimeIntegrator_->init_integrate_velocity(dt());
    ghostManager_.init_integrate_velocity(dt());
    // account for other fixes doing time integration
    interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_VELOCITY);

    atomTimeIntegrator_->init_integrate_position(dt());
    ghostManager_.init_integrate_position(dt());
    // account for other fixes doing time integration
    interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_POSITION);
  }
  //-------------------------------------------------------------------
  void ATC_Method::post_init_integrate()
  {
    ghostManager_.post_init_integrate();
  }
  //-------------------------------------------------------------------
  void ATC_Method::pre_exchange()
  {
    adjust_xref_pbc();
    // call interscale manager to sync atc per-atom data with lammps array ahead of parallel communication
    interscaleManager_.prepare_exchange();

    // change types based on moving from internal region to ghost region
    if ((atomToElementMapType_ == EULERIAN) && (step() % atomToElementMapFrequency_ == 0)) {
      ghostManager_.pre_exchange();
    }
  }
  //-------------------------------------------------------------------
  void ATC_Method::setup_pre_exchange()
  {
    adjust_xref_pbc();
    // call interscale manager to sync atc per-atom data with lammps array ahead of parallel communication
    interscaleManager_.prepare_exchange();
  }
  //-------------------------------------------------------------------
  void ATC_Method::pre_neighbor()
  {
    // reset quantities arising from atom exchange
    reset_nlocal();

    interscaleManager_.post_exchange();

    // forward_comm should go here
  }
  //-------------------------------------------------------------------
  void ATC_Method::min_post_force()
  {
    post_force();
  }
  //-------------------------------------------------------------------
  void ATC_Method::post_force()
  {
    // this resets allow for the possibility of other fixes modifying positions and velocities, e.g. walls, but reduces efficiency
    interscaleManager_.lammps_force_reset();
  }
  //--------------------------------------------------------
  void ATC_Method::final_integrate()
  {
    atomTimeIntegrator_->final_integrate(dt());
    ghostManager_.final_integrate(dt());
    // account for other fixes doing time integration
    interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_VELOCITY);
  }
  //-------------------------------------------------------------------
  void ATC_Method::post_final_integrate()
  {
    if (atomicWeightsWriteFlag_ && (step() % atomicWeightsWriteFrequency_ == 0)) {
      write_atomic_weights(atomicWeightsFile_,atomVolume_->quantity());
    }
  }
  //-------------------------------------------------------------------
  void ATC_Method::end_of_step()
  {
    localStep_ += 1;
  }
  //--------------------------------------------------------------
  void ATC_Method::finish()
  {
    // FE Engine
    if (feEngine_) feEngine_->finish();
    feEngine_->departition_mesh();
  }

  //--------------------------------------------------------------
  /** method to add new fields to the included list */
  //--------------------------------------------------------------
  void ATC_Method::add_fields(map<FieldName,int> & newFieldSizes)
  {
    map<FieldName,int>::const_iterator field;
    for (field = newFieldSizes.begin(); field!=newFieldSizes.end(); field++) {
      FieldName thisField = field->first;
      int thisSize = field->second;
      if (fieldSizes_.find(thisField)==fieldSizes_.end()) {
          fieldSizes_[thisField] = thisSize;
      }
    }
  }

//-------------------------------------------------------------------
  void ATC_Method::set_reference_potential_energy(void)
  {
    if (setRefPE_) {
      if (setRefPEvalue_) {
        nodalRefPotentialEnergy_->set_quantity() = refPEvalue_;
        setRefPEvalue_ = false;
      }
      else if (readRefPE_) {
        if (LammpsInterface::instance()->rank_zero()) {
          stringstream ss;
          ss << "reading reference potential energy from " << nodalRefPEfile_;
          LammpsInterface::instance()->print_msg(ss.str());
        }
        (nodalRefPotentialEnergy_->set_quantity()).from_file(nodalRefPEfile_);
        readRefPE_ = false;
      }
      else {
        hasRefPE_ = false;
        SPAR_MAN * referenceAccumulant = interscaleManager_.sparse_matrix("ReferenceAccumulant");
        if (referenceAccumulant) {
          referenceAccumulant->set_quantity() = accumulant_->quantity();
        }
        DIAG_MAN * referenceAccumulantInverseVolumes = interscaleManager_.diagonal_matrix("ReferenceAccumulantInverseVolumes");
        if (referenceAccumulantInverseVolumes) {
          referenceAccumulantInverseVolumes->set_quantity() = accumulantInverseVolumes_->quantity();
        }
        PAQ * atomicRefPe = interscaleManager_.per_atom_quantity("AtomicReferencePotential");
        if (!atomicRefPe) {
          throw ATC_Error("ATC_Method::set_reference_potential_energy - atomic reference PE object was not created during construct_transfers");
        }
        PAQ* pe = interscaleManager_.per_atom_quantity("AtomicPotentialEnergy");
        if (!pe) {
          throw ATC_Error("ATC_Method::set_reference_potential_energy - atomic PE object was not created during construct_transfers");
        }
        atomicRefPe->set_quantity() = pe->quantity();
        atomicRefPe->fix_quantity();
      }
      setRefPE_ = false;
      hasRefPE_ = true;
    }
  }
//-------------------------------------------------------------------


  //=================================================================
  // memory management and processor information exchange
  //=================================================================


  //-----------------------------------------------------------------
  // number of doubles
  //-----------------------------------------------------------------
  int ATC_Method::doubles_per_atom() const
  {

    int doubles = 4;
    doubles += interscaleManager_.memory_usage();
    return doubles;
  }

  //-----------------------------------------------------------------
  // memory usage of local atom-based arrays
  //-----------------------------------------------------------------
  int ATC_Method::memory_usage()
  {
    int bytes = doubles_per_atom();
    bytes *= lammpsInterface_->nmax() * sizeof(double);
    return bytes;
  }

  //-----------------------------------------------------------------
  // allocate local atom-based arrays
  //-----------------------------------------------------------------
  void ATC_Method::grow_arrays(int nmax)
  {
    xref_ =
      lammpsInterface_->grow_2d_double_array(xref_,nmax,3,"fix_atc:xref");

    perAtomOutput_ =
      lammpsInterface_->grow_2d_double_array(perAtomOutput_,nmax,sizePerAtomCols_,"fix_atc:perAtomOutput");
    interscaleManager_.grow_arrays(nmax);
  }

  //-----------------------------------------------------------------
  // copy values within local atom-based arrays
  //-----------------------------------------------------------------
  void ATC_Method::copy_arrays(int i, int j)
  {
    xref_[j][0] = xref_[i][0];
    xref_[j][1] = xref_[i][1];
    xref_[j][2] = xref_[i][2];

    for (int ii = 0 ; ii < sizePerAtomCols_ ; ii++ ) {
      perAtomOutput_[j][ii] = perAtomOutput_[i][ii];
    }
    interscaleManager_.copy_arrays(i,j);
  }

  //-----------------------------------------------------------------
  // pack values in local atom-based arrays for exchange with another proc
  //-----------------------------------------------------------------
  int ATC_Method::pack_exchange(int i, double *buf)
  {
    buf[0] = xref_[i][0];
    buf[1] = xref_[i][1];
    buf[2] = xref_[i][2];

    int j = 4;
    for (int ii = 0 ; ii < sizePerAtomCols_ ; ii++ ) {
      buf[j++] = perAtomOutput_[i][ii];
    }
    int interscaleSizeComm = interscaleManager_.pack_exchange(i,&buf[j]);
    return sizeComm_ + interscaleSizeComm;
  }

  //-----------------------------------------------------------------
  // unpack values in local atom-based arrays from exchange with another proc
  //-----------------------------------------------------------------
  int ATC_Method::unpack_exchange(int nlocal, double *buf)
  {
    xref_[nlocal][0] = buf[0];
    xref_[nlocal][1] = buf[1];
    xref_[nlocal][2] = buf[2];

    int j = 4;
    for (int ii = 0 ; ii < sizePerAtomCols_ ; ii++ ) {
      perAtomOutput_[nlocal][ii] = buf[j++];
    }
    int interscaleSizeComm = interscaleManager_.unpack_exchange(nlocal,&buf[j]);
    return sizeComm_ + interscaleSizeComm;
  }

  //-----------------------------------------------------------------
  // pack values in local atom-based arrays from exchange with another proc
  //-----------------------------------------------------------------
  int ATC_Method::pack_comm(int n, int *list, double *buf,
                            int pbc_flag, int *pbc)
  {
    int i,j,m;
    double dx = 0,dy = 0,dz = 0;

    int * num_bond = lammpsInterface_->num_bond();
    int ** bond_atom = lammpsInterface_->bond_atom();

    m = 0;
    if (pbc_flag == 0) {
      for (i = 0; i < n; i++) {
        j = list[i];
        buf[m++] = xref_[j][0];
        buf[m++] = xref_[j][1];
        buf[m++] = xref_[j][2];

        if (num_bond) {
          buf[m++] = num_bond[j];
          for (int ii = 0; ii < lammpsInterface_->bond_per_atom(); ii++) {
            buf[m++] = bond_atom[j][ii];
          }
        }
      }
    }
    else {
      if (lammpsInterface_->domain_triclinic() == 0) {
        dx = pbc[0]*Xprd_;
        dy = pbc[1]*Yprd_;
        dz = pbc[2]*Zprd_;
      }
      else {
        dx = pbc[0]*Xprd_ + pbc[5]*XY_ + pbc[4]*XZ_;
        dy = pbc[1]*Yprd_ + pbc[3]*YZ_;
        dz = pbc[2]*Zprd_;
      }
      for (i = 0; i < n; i++) {
        j = list[i];
        buf[m++] = xref_[j][0] + dx;
        buf[m++] = xref_[j][1] + dy;
        buf[m++] = xref_[j][2] + dz;

        if (num_bond) {
          buf[m++] = num_bond[j];
          for (int ii = 0; ii < lammpsInterface_->bond_per_atom(); ii++) {
            buf[m++] = bond_atom[j][ii];
          }
        }

      }
    }
    return m;  // total amount of data sent
  }

  //-----------------------------------------------------------------
  // unpack values in local atom-based arrays from exchange with another proc
  //-----------------------------------------------------------------
  void ATC_Method::unpack_comm(int n, int first, double *buf)
  {
    int i,m,last;

    int * num_bond = lammpsInterface_->num_bond();
    int ** bond_atom = lammpsInterface_->bond_atom();

    m = 0;
    last = first + n;
    for (i = first; i < last; i++) {
      xref_[i][0] = buf[m++];
      xref_[i][1] = buf[m++];
      xref_[i][2] = buf[m++];

      if (num_bond) {
        num_bond[i] = static_cast<int>(buf[m++]);
        for (int ii = 0; ii < lammpsInterface_->bond_per_atom(); ii++) {
          bond_atom[i][ii] = static_cast<int>(buf[m++]);
        }
      }

    }
  }

  //-----------------------------------------------------------------
  //
  //-----------------------------------------------------------------
  int ATC_Method::comm_forward()
  {
    int size = 3;
    if (lammpsInterface_->num_bond())
      { size += lammpsInterface_->bond_per_atom()+1; }
    return size;
  }

  //-----------------------------------------------------------------
  //
  //-----------------------------------------------------------------
  void ATC_Method::reset_nlocal()
  {
    nLocalTotal_ = lammpsInterface_->nlocal();
    const int * mask = lammpsInterface_->atom_mask();
    nLocal_ = 0;
    nLocalGhost_ = 0;

    for (int i = 0; i < nLocalTotal_; ++i) {
      if (mask[i] & groupbit_) nLocal_++;
      if (mask[i] & groupbitGhost_) nLocalGhost_++;
    }

    // set up internal & ghost maps

    if (nLocal_>0) {
      // set map
      internalToAtom_.resize(nLocal_);
      int j = 0;
      // construct internalToAtom map
      //  : internal index -> local lammps atom index
      for (int i = 0; i < nLocalTotal_; ++i) {
        if (mask[i] & groupbit_) internalToAtom_(j++) = i;
      }
#ifdef EXTENDED_ERROR_CHECKING
      stringstream ss;
      ss << "Nlocal = " << nLocal_ << " but only found " << j << "atoms";
      if (j!=nLocal_) throw ATC_Error(ss.str());
#endif
      // construct reverse map
      atomToInternal_.clear();
      for (int i = 0; i < nLocal_; ++i) {
        atomToInternal_[internalToAtom_(i)] = i;
      }
    }
    if (nLocalGhost_>0) {
      // set map
      ghostToAtom_.resize(nLocalGhost_);
      int j = 0;
      for (int i = 0; i < nLocalTotal_; ++i) {
        if (mask[i] & groupbitGhost_) ghostToAtom_(j++) = i;
      }
    }

    //WIP_JAT this should not be needed at all, but a memory problem with sparse matrices requires them to be reset (possibly related to note in SparseMatrix-inl.h::_delete())
    interscaleManager_.reset_nlocal();

  }

  //-------------------------------------------------------------------
  void ATC_Method::reset_coordinates()
  {
    // update coarse graining positions for internal and ghost atoms
    atomCoarseGrainingPositions_->unfix_quantity();
    atomCoarseGrainingPositions_->quantity();
    atomCoarseGrainingPositions_->fix_quantity();
    if (atomGhostCoarseGrainingPositions_) {
      atomGhostCoarseGrainingPositions_->unfix_quantity();
      atomGhostCoarseGrainingPositions_->quantity();
      atomGhostCoarseGrainingPositions_->fix_quantity();
    }
     if (atomProcGhostCoarseGrainingPositions_) {
      atomProcGhostCoarseGrainingPositions_->unfix_quantity();
      atomProcGhostCoarseGrainingPositions_->quantity();
      atomProcGhostCoarseGrainingPositions_->fix_quantity();
    }
  }

  //-----------------------------------------------------------------
  //
  //-----------------------------------------------------------------
  void ATC_Method::write_atomic_weights(const string filename, const DIAG_MAT & atomicVolumeMatrix)
  {
    int nlocal = lammpsInterface_->nlocal();
    int nlocalmax;
    LammpsInterface::instance()->int_allmax(&nlocal,&nlocalmax);
    int natoms = int(lammpsInterface_->natoms());
    ofstream out;
    const char* fname = &filename[0];

    // create tag to local id map for this processor
    map <int,int> id2tag;
    map <int,int>::const_iterator itr;
    int * atom_tag = lammpsInterface_->atom_tag();
    for (int i = 0; i < nlocal; ++i) {
      id2tag[i] = atom_tag[i];
    }

    int comm_rank = LammpsInterface::instance()->comm_rank();
    int nprocs;
    LammpsInterface::instance()->int_allmax(&comm_rank,&nprocs);
    nprocs += 1;

    if (comm_rank == 0) {
      out.open(fname);
      // print header lines
      out << "Atomic Weights for LAMMPS/atc analysis\n";
      out << " \n";
      out << natoms << " Atoms in system\n";
      out << " \n";
      // print atomic weights from proc 0
      for(int i = 0; i < nlocal; i++) {
        out << id2tag[i] << "  " << atomicVolumeMatrix(i,i) << "\n";
      }
    }

    if (nprocs > 1) {
      int max_size,send_size;
      send_size = nlocal;
      LammpsInterface::instance()->int_allmax(&send_size,&max_size);

      if (comm_rank == 0) {
        int *intbuf = new int[max_size];
        double *buf = new double[max_size];
        for (int iproc = 1; iproc < nprocs; iproc++) {
          LammpsInterface::instance()->int_recv(intbuf,max_size,iproc);
          LammpsInterface::instance()->recv(buf,max_size,iproc);
          for (int i = 0; i < max_size; i++) {
            out << intbuf[i] << "  " << buf[i] << "\n";
          }
        }
        delete[] intbuf;
        delete[] buf;
      } else {
        int *intbuf = new int[send_size];
        double *buf = new double[send_size];
        for (int i = 0; i < send_size; i++) {
          intbuf[i] = id2tag[i];
          buf[i] = atomicVolumeMatrix(i,i);
        }
        LammpsInterface::instance()->int_send(intbuf,send_size);
        LammpsInterface::instance()->send(buf,send_size);
        delete[] intbuf;
        delete[] buf;
      }
    }

    if (comm_rank == 0) {
      out.close();
    }
  }

  //-----------------------------------------------------------------
  //
  //-----------------------------------------------------------------
  void ATC_Method::compute_consistent_md_mass_matrix(const SPAR_MAT & shapeFunctionMatrix,
                                                       SPAR_MAT & mdMassMatrix) const
  {

    int nCols = shapeFunctionMatrix.nCols();
    DENS_MAT massMatrixLocal(nCols,nCols);
    DENS_MAT denseMdMassMatrix(nCols,nCols);
    if (nLocal_>0)
      massMatrixLocal = shapeFunctionMatrix.transMat(shapeFunctionMatrix);

    lammpsInterface_->allsum(massMatrixLocal.ptr(),
                             denseMdMassMatrix.ptr(),
                             denseMdMassMatrix.size());
    mdMassMatrix.reset(denseMdMassMatrix,1.e-10);
  }

  //=================================================================
  // Interscale operators
  //=================================================================
  // in the spirit of the current design of ATC: atoms local, nodes global




  bool ATC_Method::nodal_influence(const int groupbit,
                              set<int> & nset, set<int> & aset, double tol)
  {
    int nghost = nodal_influence(groupbit,nset,aset,true,tol);
    int local_nghost = nghost;
    lammpsInterface_->int_allsum(&local_nghost,&nghost);
    if (nghost == 0) {
       nodal_influence(groupbit,nset,aset,false,tol);
    }
    return (nghost > 0);
  }
  int ATC_Method::nodal_influence(const int groupbit,
        set<int> & nset, set<int> & aset, bool ghost, double tol)
  {
    Array<int> & amap = (ghost) ? ghostToAtom_ : internalToAtom_;
    int natoms = (ghost) ? nLocalGhost_ : nLocal_;
    DENS_MAT influence(nNodes_,1);
    DENS_MAT atomInfluence(natoms,1);
    const int *mask = lammpsInterface_->atom_mask();
    for (int i = 0; i < natoms; i++) {
      if (mask[amap(i)] & groupbit) {
         atomInfluence(i,0) = 1;
         aset.insert(i);
      }
    }
    // relies on shape functions
    if (ghost) {
      restrict_volumetric_quantity(atomInfluence,influence,(interscaleManager_.per_atom_sparse_matrix("InterpolantGhost"))->quantity());
    }
    else {
    restrict_volumetric_quantity(atomInfluence,influence);
    }

    DENS_MAT localInfluence = influence;
    lammpsInterface_->allsum(localInfluence.ptr(),
                             influence.ptr(),
                             influence.size());

    for (int i = 0; i < nNodes_; i++) {
      if (fabs(influence(i,0)) > tol)  { nset.insert(i);  }
    }
    return aset.size();
  }


  //--------------------------------------------------------
  void ATC_Method::restrict_volumetric_quantity(const MATRIX & atomData,
                                                MATRIX & nodeData,
                                                const SPAR_MAT & shpFcn)
  {
    // computes nodeData = N*DeltaVAtom*atomData where N are the shape functions
    DENS_MAT workNodeArray(nodeData.nRows(),nodeData.nCols());
    //DENS_MAT workNodeArray;


    if (atomData.nRows() > 0) { // or shpFcn_???
      workNodeArray = shpFcn.transMat(atomData);
    }
    int count = nodeData.nRows()*nodeData.nCols();
    lammpsInterface_->allsum(workNodeArray.ptr(),nodeData.ptr(),count);
    return;
  }


  //--------------------------------------------------------
  void ATC_Method::restrict_volumetric_quantity(const MATRIX & atomData,
                                                MATRIX & nodeData)
  {
    restrict_volumetric_quantity(atomData,nodeData,shpFcn_->quantity());
    return;
  }


  //--------------------------------------------------------
  void ATC_Method::prolong(const MATRIX & nodeData,
                             MATRIX & atomData)
  {
    // computes the finite element interpolation at atoms atomData = N*nodeData
    if (nLocal_>0) {
      atomData = (shpFcn_->quantity())*nodeData;
    }
    return;
  }

  //========================================================
  // FE functions
  //========================================================

  //--------------------------------------------------------
  void ATC_Method::output()
  {
//  if (lammpsInterface_->comm_rank() == 0) {
      compute_nodeset_output();
      compute_faceset_output();
      compute_elementset_output();
//  }
  }
  //--------------------------------------------------------
  void ATC_Method::compute_nodeset_output(void)
  {
    map< pair <string, FieldName>, NodesetOperationType >::const_iterator iter;
    for (iter = nsetData_.begin(); iter != nsetData_.end();iter++){
      pair <string, FieldName> id = iter->first;
      string nsetName = id.first;
      FieldName field = id.second;
      double sum = 0.0;
      const set<int> nset = feEngine_->fe_mesh()->nodeset(nsetName);
      const DENS_MAT & thisField = (fields_.find(field)->second).quantity();
      set< int >::const_iterator itr;
      for (itr = nset.begin(); itr != nset.end();itr++){
        int node = *itr;
        sum += thisField(node,0);
      }
      string name = nsetName + "_" + field_to_string(field);
      if (iter->second == NODESET_AVERAGE) {
        sum /= nset.size();
        name = "average_"+name;
      }
      feEngine_->add_global(name, sum);
    }
  }
  //--------------------------------------------------------
  void ATC_Method::compute_faceset_output(void)
  {
    map < pair<string,string>, FacesetIntegralType >::const_iterator iter;
    DENS_MAT values;
    for (iter = fsetData_.begin(); iter !=  fsetData_.end(); iter++) {
      string bndyName  = (iter->first).first;
      string fieldName = (iter->first).second;
      const set< PAIR > & faceSet = (feEngine_->fe_mesh())->faceset(bndyName);
      ATOMIC_DATA::iterator data_iter = filteredData_.find(fieldName);
      if (data_iter == filteredData_.end()) {
        string msg = "Specified fieldName "+fieldName+
          " not found in filteredData_ while attempting surface integration";
        throw ATC_Error(msg);
      }
      const DENS_MAT & data =  ((data_iter->second).quantity());
      string stem = bndyName + "_" + fieldName + "_";
      bool tf = false;
      if (iter->second == CONTOUR_INTEGRAL) {
        stem = "contour_"+stem;
        tf = true;
      }
      feEngine_->field_surface_flux(data,faceSet,values,tf);
      for (int i = 0; i < values.nRows() ; ++i ) {
        string name = stem + to_string(i+1);
        feEngine_->add_global(name, values(i,0));
      }
    }
  }
  //--------------------------------------------------------
  void ATC_Method::compute_elementset_output(void)
  {
    map< pair <string, FieldName>, ElementsetOperationType >::const_iterator iter;
    for (iter = esetData_.begin(); iter != esetData_.end();iter++){
      pair <string, FieldName> id = iter->first;
      string esetName = id.first;
      FieldName field = id.second;
      const ESET eset = feEngine_->fe_mesh()->elementset(esetName);
      const DENS_MAT & thisField = (fields_.find(field)->second).quantity();
      DENS_VEC total = feEngine_->integrate(thisField,eset);
      string name = esetName + "_" + field_to_string(field);
      if (iter->second == ELEMENTSET_AVERAGE) {
        DENS_MAT ones(nNodes_,0); ones = 1;
        DENS_VEC V = feEngine_->integrate(ones,eset);
        total /= V[0];
        name = "average_"+name;
      }
      if (total.size() == 1) {
        feEngine_->add_global(name, total[0]);
      }
      else {
        for (int i = 0; i < total.size(); i++) {
          feEngine_->add_global(name+to_string(i), total[i]);
        }
      }
    }
  }


  //=================================================================
  //
  //=================================================================
  //--------------------------------------------------------
  bool ATC_Method::read_atomic_ref_positions(const char * filename)
  {
    int nlocal = lammpsInterface_->nlocal();
    ifstream in;
    const int lineSize = 256;
    char line[lineSize];

    // create tag to local id map for this processor
    map <int,int> tag2id;
    map <int,int>::const_iterator itr;
    int * atom_tag = lammpsInterface_->atom_tag();
    for (int i = 0; i < nlocal; ++i) {
      tag2id[atom_tag[i]] = i;
    }

    // get number of atoms
    int natoms = 0;
    int ncols = 0;
    int style = LammpsInterface::CHARGE_STYLE;
    if (LammpsInterface::instance()->rank_zero()) {
      in.open(filename);
      string msg;
      string name = filename;
      msg = "no "+name+" reference file found";
      if (! in.good()) throw ATC_Error(msg);
      in.getline(line,lineSize); // header
      in.getline(line,lineSize); // blank line
      in >> natoms;
      stringstream ss;
      ss << "found " << natoms << " atoms in reference " << filename ;
      while(in.good()) {
        in.getline(line,lineSize);
        string str(line);
        int pos = str.find("Atoms");
        if (pos > -1) {
          in.getline(line,lineSize); // blank line
          break;
        }
      }
      in.getline(line,lineSize);
      std::vector<std::string> tokens;
      ATC_Utility::command_strings(line, tokens);
      ncols = tokens.size();
      switch (ncols) {
      // atomic: id type x y z
      case 5:
      case 8:
         ss << " style:atomic";
         style = LammpsInterface::ATOMIC_STYLE;
         break;
      // charge: id type q x y z
      // molecule : id molecule-ID type x y z
      case 6:
         ss << " style:charge";
         style = LammpsInterface::CHARGE_STYLE;
         break;
      // full  : id molecule-ID type q x y z
      case 7:
         ss << " style:full";
         style = LammpsInterface::FULL_STYLE;
         break;
      default:
         throw ATC_Error("cannot determine atom style, columns:"+to_string(ncols));
         break;
      }
      LammpsInterface::instance()->print_msg(ss.str());
      in.close();
    }
    LammpsInterface::instance()->int_broadcast(&natoms);

    // read atoms and assign
    if (LammpsInterface::instance()->rank_zero()) {
      in.open(filename);
      while(in.good()) {
        in.getline(line,lineSize);
        string str(line);
        int pos = str.find("Atoms");
        if (pos > -1) {
          in.getline(line,lineSize); // blank line
          break;
        }
      }
    }
    int nread = 0,type = -1, tag = -1, count = 0, mId = -1;
    double x[3]={0,0,0}, q =0;
    while (nread < natoms) {
      if (LammpsInterface::instance()->rank_zero()) {
         in.getline(line,lineSize);
         stringstream ss (line,stringstream::in | stringstream::out);
         if      (style == LammpsInterface::CHARGE_STYLE)
         ss >> tag >> type >> q >> x[0] >> x[1] >> x[2];
         else if (style == LammpsInterface::FULL_STYLE)
         ss >> tag >> mId >> type >> q >> x[0] >> x[1] >> x[2];
         else
         ss >> tag >> type >> x[0] >> x[1] >> x[2];
         nread++;
      }
      LammpsInterface::instance()->int_broadcast(&nread);
      LammpsInterface::instance()->int_broadcast(&tag);
      LammpsInterface::instance()->broadcast(x,3);
      itr = tag2id.find(tag);
      if (itr != tag2id.end()) {
        int iatom = itr->second;
        xref_[iatom][0] = x[0];
        xref_[iatom][1] = x[1];
        xref_[iatom][2] = x[2];
        count++;
      }
    }
    if (LammpsInterface::instance()->rank_zero()) {
      in.close();
      stringstream ss;
      ss << "read  " << natoms << " reference positions";
      LammpsInterface::instance()->print_msg(ss.str());
    }
    if (count != nlocal)
       throw ATC_Error("reset "+to_string(count)+" atoms vs "+to_string(nlocal));


    return true;
  }

//--------------------------------------------------------
  void ATC_Method::remap_ghost_ref_positions(void)
  {

    int nlocal = lammpsInterface_->nlocal();
    int nghost = lammpsInterface_->nghost();

    double box_bounds[2][3];
    lammpsInterface_->box_bounds(box_bounds[0][0],box_bounds[1][0],
                                    box_bounds[0][1],box_bounds[1][1],
                                    box_bounds[0][2],box_bounds[1][2]);
    double xlo = box_bounds[0][0], xhi = box_bounds[1][0];
    double ylo = box_bounds[0][1], yhi = box_bounds[1][1];
    double zlo = box_bounds[0][2], zhi = box_bounds[1][2];

    double box_length[3];
    for (int k = 0; k < 3; k++) {
      box_length[k] = box_bounds[1][k] - box_bounds[0][k];
    }
    double Lx = box_length[0], Ly = box_length[1], Lz = box_length[2];

    // create tag to local id map for this processor
    map <int,int> tag2id;
    map <int,int>::const_iterator itr;
    int * atom_tag = lammpsInterface_->atom_tag();
    for (int i = 0; i < nlocal; ++i) {
      tag2id[atom_tag[i]] = i;
    }

    // loop over ghosts
    double ** x = lammpsInterface_->xatom();
    for (int j = nlocal; j < nlocal + nghost; j++) {
      int tag = atom_tag[j];
      int i = tag2id[tag];
      //itr = tag2id.find(tag);
      //if (itr != tag2id.end())
      double* xj = x[j];
      double* Xj = xref_[j];
      //double Xj[3];
      double* Xi = xref_[i];
      // the assumption is that xref_[j] has been shuffled
      // so make an image of xref_[i] that is close to x[j]
      if (xj[0] <= xlo) Xj[0] = Xi[0] -Lx;
      if (xj[0] >= xhi) Xj[0] = Xi[0] +Lx;
      if (xj[1] <= ylo) Xj[1] = Xi[1] -Ly;
      if (xj[1] >= yhi) Xj[1] = Xi[1] +Ly;
      if (xj[2] <= zlo) Xj[2] = Xi[2] -Lz;
      if (xj[2] >= zhi) Xj[2] = Xi[2] +Lz;
    }
  }
};