// ATC transfer headers
#include "PerAtomQuantityLibrary.h"
#include "ATC_Transfer.h"
#include "FE_Engine.h"
#include "LammpsInterface.h"
#include <typeinfo>
#include <sstream>
#include <iostream>

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

using ATC_Utility::to_string;

namespace ATC {

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomToElementMap
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomToElementMap::AtomToElementMap(ATC_Method * atc,
                                     PerAtomQuantity<double> * atomPositions,
                                     AtomType atomType) :
    ProtectedAtomQuantity<int>(atc,1,atomType),
    atomPositions_(atomPositions)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomPositions_)
      atomPositions_ = interscaleManager.per_atom_quantity("AtomicCoarseGrainingPositions");
    if (!atomPositions_)
      throw ATC_Error("AtomToElementMap::AtomToElementMap - atom position quantity is undefined");

    atomPositions_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomToElementMap::~AtomToElementMap()
  {
    atomPositions_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomToElementMap::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<int>::reset();
      const DENS_MAT & position(atomPositions_->quantity());
      const FE_Mesh * feMesh = atc_.fe_engine()->fe_mesh();
      int nsd = atc_.nsd();
      DENS_VEC coords(nsd);
      for (int i = 0; i < quantity_.nRows(); i++) {
        for (int j = 0; j < nsd; j++) {
          coords(j) = position(i,j);
        }
        quantity_(i,0) = feMesh->map_to_element(coords);
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomInElementSet
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomInElementSet::AtomInElementSet(ATC_Method * atc,
                     PerAtomQuantity<int> * map,
                     ESET eset, int type):
  atc_(atc,INTERNAL), 
  map_(map),eset_(eset),type_(type),
  quantityToLammps_(atc_.atc_to_lammps_map())
  {
    map_->register_dependence(this);
    needReset_ = true;
    
  }
  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomInElementSet::~AtomInElementSet()
  {
    map_->remove_dependence(this);
  }
  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomInElementSet::reset() 
  {
    if (map_->need_reset() || needReset_) {
      list_.clear();
      INT_ARRAY map = map_->quantity();
      int * type = ATC::LammpsInterface::instance()->atom_type(); 
      ESET::const_iterator itr;
      const Array<int> & quantityToLammps = atc_.atc_to_lammps_map();
      for (int i = 0; i < map.size(); i++) {
        for (itr=eset_.begin(); itr != eset_.end(); itr++) {
          if (map(i,0) == *itr) {
            int a = quantityToLammps(i);
            if (type[a] == type_) {
              list_.push_back(ID_PAIR(i,a)); 
              break;
            }
          }
        }
      }
      needReset_ = false;
    }
  }
  //--------------------------------------------------------
  //  qauntity
  //--------------------------------------------------------
  const ID_LIST &  AtomInElementSet::quantity()
  {
    reset();
    return list_;
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomVolumeUser
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomVolumeUser::AtomVolumeUser(ATC_Method * atc,
                                 map<int,double> & atomGroupVolume,
                                 AtomType atomType) :
    ProtectedAtomDiagonalMatrix<double>(atc,atomType),
    atomGroupVolume_(atomGroupVolume),
    lammpsInterface_(atc->lammps_interface()),
    atcToLammps_(atc->internal_to_atom_map())
  {
    // do nothing
  }
  
  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomVolumeUser::reset() const
  {
    if (need_reset()) {
      PerAtomDiagonalMatrix<double>::reset();
      const int *mask = lammpsInterface_->atom_mask();
      quantity_ = 0.;

      map<int, double>::const_iterator igroup;
      for (igroup = atomGroupVolume_.begin(); igroup != atomGroupVolume_.end(); igroup++) {
        int gid = igroup->first;
        double weight = igroup->second;
        for (int i = 0; i < quantity_.nRows(); ++i) {
          if (mask[atcToLammps_(i)] & gid) {
            quantity_(i,i) = weight;
          }
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomVolumeGroup
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomVolumeGroup::AtomVolumeGroup(ATC_Method * atc,
                                   map<int,double> & atomGroupVolume,
                                   AtomType atomType) :
    AtomVolumeUser(atc,atomGroupVolume,atomType),
    atomGroupVolume_(atomGroupVolume),
    lammpsInterface_(atc->lammps_interface()),
    atcToLammps_(atc->internal_to_atom_map())
  {
    // Uses group bounding box as total volume
    // ASSUME will *only* work if atoms exactly fill the box
    int ngroup = lammpsInterface_->ngroup();
    const int *mask = lammpsInterface_->atom_mask();
    double * bounds;
      
    bounds = new double[6];
    for (int i = 0; i < ngroup; ++i) {
      lammpsInterface_->group_bounds(i, bounds);
      atomGroupVolume_[lammpsInterface_->group_bit(i)] = 
        (bounds[1]-bounds[0])*(bounds[3]-bounds[2])*(bounds[5]-bounds[4]);
    }
    delete [] bounds;
      
    INT_VECTOR localCount(ngroup);
    INT_VECTOR globalCount(ngroup);
      
    // loop over atoms
    localCount = 0;
    for (int i = 0; i < atcToLammps_.size(); ++i) {
      for (int j = 0; j < ngroup; ++j) {
        if (mask[atcToLammps_(i)] & lammpsInterface_->group_bit(j))
          localCount(j)++;
      }
    }
      
    // communication to get total atom counts per group
    lammpsInterface_->int_allsum(localCount.ptr(),
                                 globalCount.ptr(),ngroup);
    
    for (int i = 0; i < ngroup; ++i) {
      int iGroupBit = lammpsInterface_->group_bit(i); 
      if (globalCount(i) > 0)
        atomGroupVolume_[iGroupBit] /= globalCount(i);
      else
        atomGroupVolume_[iGroupBit] = 0;
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomVolumeLattice
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomVolumeLattice::AtomVolumeLattice(ATC_Method * atc,
                                       AtomType atomType) :
    ProtectedAtomDiagonalMatrix<double>(atc,atomType),
    lammpsInterface_(atc->lammps_interface())
  {
    // do nothing
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomVolumeLattice::reset() const
  {
    if (need_reset()) {
      PerAtomDiagonalMatrix<double>::reset();
      quantity_ = lammpsInterface_->volume_per_atom();
    }
  }
  
  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomVolumeElement
  //--------------------------------------------------------
  //--------------------------------------------------------
  
  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomVolumeElement::AtomVolumeElement(ATC_Method * atc,
                                       PerAtomQuantity<int> * atomElement,
                                       AtomType atomType) :
    ProtectedAtomDiagonalMatrix<double>(atc,atomType),
    atomElement_(atomElement),
    lammpsInterface_(atc->lammps_interface()),
    feMesh_((atc->fe_engine())->fe_mesh())
  {
    if (!atomElement_) {
      InterscaleManager & interscaleManager(atc->interscale_manager());
      atomElement_ = interscaleManager.per_atom_int_quantity("AtomElement");
    }
    atomElement_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  Destructor
  //--------------------------------------------------------
  AtomVolumeElement::~AtomVolumeElement()
  {
    atomElement_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomVolumeElement::reset() const
  {
    if (need_reset()) {
      PerAtomDiagonalMatrix<double>::reset();
      int nElts = feMesh_->num_elements();
      int nLocal = quantity_.nRows();
      _elementAtomCountLocal_.reset(nElts);
      _elementAtomCount_.resize(nElts);
      const INT_ARRAY & atomElement(atomElement_->quantity());
       _elementAtomVolume_.resize(nElts);
      
      // determine number of atoms in each element, partial sum
      for (int i = 0; i < nLocal; ++i) {
        _elementAtomCountLocal_(atomElement(i,0)) += 1;
      }
      
      // mpi to determine total atoms
      lammpsInterface_->int_allsum(_elementAtomCountLocal_.ptr(),_elementAtomCount_.ptr(),nElts);
      
      // divide element volume by total atoms to get atomic volume
      if (nLocal>0) {
        for (int i = 0; i < nElts; ++i) {
          
          double minx, maxx, miny, maxy, minz, maxz;
          feMesh_->element_coordinates(i,_nodalCoordinates_);
          minx = _nodalCoordinates_(0,0); maxx = _nodalCoordinates_(0,0);
          miny = _nodalCoordinates_(1,0); maxy = _nodalCoordinates_(1,0);
          minz = _nodalCoordinates_(2,0); maxz = _nodalCoordinates_(2,0);
          for (int j = 1; j < _nodalCoordinates_.nCols(); ++j) {
            if (_nodalCoordinates_(0,j)<minx) minx = _nodalCoordinates_(0,j);
            if (_nodalCoordinates_(0,j)>maxx) maxx = _nodalCoordinates_(0,j);
            if (_nodalCoordinates_(1,j)<miny) miny = _nodalCoordinates_(1,j);
            if (_nodalCoordinates_(1,j)>maxy) maxy = _nodalCoordinates_(1,j);
            if (_nodalCoordinates_(2,j)<minz) minz = _nodalCoordinates_(2,j);
            if (_nodalCoordinates_(2,j)>maxz) maxz = _nodalCoordinates_(2,j);
          }
          double eltVol = (maxx-minx)*(maxy-miny)*(maxz-minz);
          if (eltVol<0) eltVol *= -1.;
          _elementAtomVolume_(i) = eltVol/_elementAtomCount_(i);
        }
        
        for (int i = 0; i < nLocal; ++i)
          quantity_(i,i) = _elementAtomVolume_(atomElement(i,0));
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomVolumeRegion
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomVolumeRegion::AtomVolumeRegion(ATC_Method * atc,
                                     DENS_MAN * atomCoarseGrainingPositions,
                                     AtomType atomType) :
    ProtectedAtomDiagonalMatrix<double>(atc,atomType),
    lammpsInterface_(atc->lammps_interface())
  {
    if (!atomCoarseGrainingPositions) {
      InterscaleManager & interscaleManager(atc->interscale_manager());
      atomCoarseGrainingPositions_ = interscaleManager.per_atom_quantity("AtomicCoarseGrainingPositions");
    }
    atomCoarseGrainingPositions_->register_dependence(this);

    // compute volumes and atom counts in each region
    int nregion = lammpsInterface_->nregion();
    regionalAtomVolume_.resize(nregion);
    
    for (int i = 0; i < nregion; ++i) {
      regionalAtomVolume_(i) = 
        (lammpsInterface_->region_xhi(i)-lammpsInterface_->region_xlo(i))
        *(lammpsInterface_->region_yhi(i)-lammpsInterface_->region_ylo(i))
        *(lammpsInterface_->region_zhi(i)-lammpsInterface_->region_zlo(i));
    }
      
    INT_VECTOR localCount(nregion);
    INT_VECTOR globalCount(nregion);  
    // loop over atoms
    localCount = 0;
    const DENS_MAT atomicCoordinates(atomCoarseGrainingPositions->quantity());
    for (int i = 0; i < quantity_.nRows(); ++i) {
      for (int j = 0; j < nregion; ++j) {
        if (lammpsInterface_->region_match(j,
                                           atomicCoordinates(i,0),
                                           atomicCoordinates(i,1),
                                           atomicCoordinates(i,2))) {
          localCount(j)++;
        }
      }
    }
    // communication to get total atom counts per group
    lammpsInterface_->int_allsum(localCount.ptr(),
                                 globalCount.ptr(),nregion);
        
    for (int i = 0; i < nregion; ++i) {
      if (globalCount(i) > 0)
        regionalAtomVolume_(i) /= globalCount(i);
      else
        regionalAtomVolume_(i) = 0;
    }
  }

  //--------------------------------------------------------
  //  Destructor
  //--------------------------------------------------------
  AtomVolumeRegion::~AtomVolumeRegion()
  {
    atomCoarseGrainingPositions_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomVolumeRegion::reset() const
  {
    if (need_reset()) {
      PerAtomDiagonalMatrix<double>::reset();
      const DENS_MAT & atomicCoordinates(atomCoarseGrainingPositions_->quantity());  
      for (int i = 0; i < quantity_.nRows(); i++) {
        for (int iregion = 0; iregion < lammpsInterface_->nregion(); iregion++) {
          if (lammpsInterface_->region_match(iregion,
                                             atomicCoordinates(i,0),
                                             atomicCoordinates(i,0),
                                             atomicCoordinates(i,2)))
            quantity_(i,i) = regionalAtomVolume_(iregion);
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomVolumeFile
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomVolumeFile::AtomVolumeFile(ATC_Method * atc,
                                 const string & atomVolumeFile,
                                 AtomType atomType) :
    ProtectedAtomDiagonalMatrix<double>(atc,atomType),
    atomVolumeFile_(atomVolumeFile),
    lammpsInterface_(atc->lammps_interface())
  {
    // do nothing
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomVolumeFile::reset() const
  {
    if (need_reset()) {
      PerAtomDiagonalMatrix<double>::reset();
      int nlocal = lammpsInterface_->nlocal();
      ifstream in;
      const int lineSize = 256;
      char line[lineSize];
      const char* fname = &atomVolumeFile_[0]; 

      // create tag to local id map for this processor
      map <int,int> tag2id;
      map <int,int>::const_iterator itr;
      const 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;
      if (ATC::LammpsInterface::instance()->rank_zero()) {
        in.open(fname);
        string msg;
        string name = atomVolumeFile_;
        msg = "no "+name+" atomic weights file found";
        if (! in.good()) throw ATC_Error(msg);
        in.getline(line,lineSize); // header
        in.getline(line,lineSize); // blank line
        in.getline(line,lineSize); // number of atoms 
        stringstream inss (line,stringstream::in | stringstream::out);
        inss >> natoms; // number of atoms
        stringstream ss; 
        ss << " found " << natoms << " atoms in atomic weights file";
        lammpsInterface_->print_msg(ss.str());
        if (natoms != lammpsInterface_->natoms()) { 
          throw ATC_Error("Incorrect number of atomic weights read-in!");
        }
        in.getline(line,lineSize); // blank line
      }

      // read and assign atomic weights
      int nread = 0, tag = -1, count = 0;
      double atomic_weight;
      while (nread < natoms) {
        if (ATC::LammpsInterface::instance()->rank_zero()) {
          in.getline(line,lineSize);
          stringstream ss (line,stringstream::in | stringstream::out);
          ss >> tag >> atomic_weight; 
          nread++;
        }
        lammpsInterface_->int_broadcast(&nread);
        lammpsInterface_->int_broadcast(&tag);
        lammpsInterface_->broadcast(&atomic_weight);
        itr = tag2id.find(tag);
        if (itr != tag2id.end()) {
          int iatom = itr->second;
          quantity_(iatom,0) = atomic_weight;
          count++;
        }
      }
      if (lammpsInterface_->rank_zero()) {
        in.close();
        stringstream ss; 
        ss << " read  " << nread << " atomic weights";
        lammpsInterface_->print_msg(ss.str());
      }
      if (count != nlocal) 
        throw ATC_Error("reset "+to_string(count)+" atoms vs "+to_string(nlocal));
    }
  }

  // need to add capability to take in group bit (JAT, 04/02/11)
  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomNumber
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomNumber::AtomNumber(ATC_Method * atc, AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,1,atomType),
    atc_(atc)  
  {
  }

  void AtomNumber::reset() const
  {
    int nlocal = atc_->nlocal();
    quantity_.reset(nlocal,1);
    quantity_ = 1;
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomTypeVector
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomTypeVector::AtomTypeVector(ATC_Method * atc, 
    vector<int> typeList, AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,
    typeList.size(),
    atomType),
    atc_(atc), ntypes_(ATC::LammpsInterface::instance()->ntypes()),
    typeList_(typeList)  
  {
    if (typeList_.size() == 0) throw ATC_Error("type list is empty");
    index_.resize(ntypes_,-1); 
    for (unsigned int i = 0; i < typeList_.size(); i++) {
      index_[typeList_[i]] = i;
    }
  }
  AtomTypeVector::AtomTypeVector(ATC_Method * atc, 
    vector<int> typeList, vector<int> groupList, AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,
    typeList.size()+groupList.size(),
    atomType),
    atc_(atc), ntypes_(ATC::LammpsInterface::instance()->ntypes()),
    typeList_(typeList),
    groupList_(groupList)  
  {
    if ((typeList_.size() == 0) && (groupList_.size() == 0)) throw ATC_Error("type/group lists are empty");
    // reverse map
    index_.resize(ntypes_,-1); 
    for (unsigned int i = 0; i < typeList_.size(); i++) {
      index_[typeList_[i]-1] = i;
    }
  }

  void AtomTypeVector::reset() const
  {
    if (need_reset()) {
      //PerAtomQuantity<double>::reset();
      int nlocal = atc_->nlocal();
      quantity_.reset(nlocal,typeList_.size()+groupList_.size());
      const Array<int> & quantityToLammps = (PerAtomQuantity<double>::atc_).atc_to_lammps_map();
      if (typeList_.size()) { 
        int * type = ATC::LammpsInterface::instance()->atom_type(); 
        for (int i = 0; i < nlocal; i++) {
          int a = quantityToLammps(i); 
          int index = index_[type[a]-1];
          if (index > -1) quantity_(i,index) = 1;
        }
      }
      int index = typeList_.size();
      if (groupList_.size()) { 
        const int * mask = ATC::LammpsInterface::instance()->atom_mask();
        for (unsigned int j = 0; j < groupList_.size(); j++) {
          int group = groupList_[j];
          for (int i = 0; i < nlocal; i++) {
            int a = quantityToLammps(i);
            if (mask[a] & group) quantity_(i,index) = 1;
          }
          index++;
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class XrefWrapper
  //--------------------------------------------------------
  //--------------------------------------------------------
  
  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  XrefWrapper::XrefWrapper(ATC_Method * atc, AtomType atomType) :
    ProtectedClonedAtomQuantity<double>(atc,atc->nsd(),atomType),
    atc_(atc)
  {
    // do nothing
  }

  //--------------------------------------------------------
  //  lammps_vector
  //--------------------------------------------------------
  double ** XrefWrapper::lammps_vector() const
  {
    return atc_->xref();
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomicMassWeightedDisplacement
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomicMassWeightedDisplacement::AtomicMassWeightedDisplacement(ATC_Method * atc,
    PerAtomQuantity<double> * atomPositions,
    PerAtomQuantity<double> * atomMasses,
    PerAtomQuantity<double> * atomReferencePositions,
    AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,atc->nsd(),atomType),
    atomPositions_(atomPositions),
    atomMasses_(atomMasses),
    atomReferencePositions_(atomReferencePositions)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomPositions_)
      atomPositions_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION,
                                                                       atomType);
    if (!atomMasses_)
      atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                                    atomType);
    if (!atomReferencePositions_)
      atomReferencePositions_ = interscaleManager.per_atom_quantity("AtomicCoarseGrainingPositions");

    atomPositions_->register_dependence(this);
    atomMasses_->register_dependence(this);
    atomReferencePositions_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomicMassWeightedDisplacement::~AtomicMassWeightedDisplacement()
  {
    atomPositions_->remove_dependence(this);
    atomMasses_->remove_dependence(this);
    atomReferencePositions_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomicMassWeightedDisplacement::reset() const
  {
    if (need_reset()) { 
      PerAtomQuantity<double>::reset();
      const DENS_MAT & position(atomPositions_->quantity());
      const DENS_MAT & mass(atomMasses_->quantity());
      const DENS_MAT & refPosition(atomReferencePositions_->quantity());

      // q = m * (x - xref)
      quantity_ = position;
      quantity_ -= refPosition;
      quantity_ *= mass;
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomicMomentum
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomicMomentum::AtomicMomentum(ATC_Method * atc,
                                 PerAtomQuantity<double> * atomVelocities,
                                 PerAtomQuantity<double> * atomMasses,
                                 AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,atc->nsd(),atomType),
    atomVelocities_(atomVelocities),
    atomMasses_(atomMasses)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomVelocities_)
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY,
                                                                    atomType);
    if (!atomMasses_)
      atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                                    atomType);

    atomVelocities_->register_dependence(this);
    atomMasses_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomicMomentum::~AtomicMomentum()
  {
    atomVelocities_->remove_dependence(this);
    atomMasses_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomicMomentum::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & mass(atomMasses_->quantity());
      const DENS_MAT & velocity(atomVelocities_->quantity());

      // q = m * v
      quantity_ = velocity;
      quantity_ *= mass;    
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class TwiceKineticEnergy
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  TwiceKineticEnergy::TwiceKineticEnergy(ATC_Method * atc,
                                         PerAtomQuantity<double> * atomVelocities,
                                         PerAtomQuantity<double> * atomMasses,
                                         AtomType atomType) :
    AtomicEnergyForTemperature(atc,atomType),
    atomVelocities_(atomVelocities),
    atomMasses_(atomMasses)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomVelocities_)
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY,
                                                                        atomType);
    if (!atomMasses_)
      atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                                    atomType);

    atomVelocities_->register_dependence(this);
    atomMasses_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  TwiceKineticEnergy::~TwiceKineticEnergy()
  {
    atomVelocities_->remove_dependence(this);
    atomMasses_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void TwiceKineticEnergy::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & mass(atomMasses_->quantity());
      const DENS_MAT & velocity(atomVelocities_->quantity());
      
      // q = m * (v dot v)
      for (int i = 0; i < quantity_.nRows(); i++) {
          quantity_(i,0) = velocity(i,0)*velocity(i,0);
        for (int j = 1; j < velocity.nCols(); j++) {
          quantity_(i,0) += velocity(i,j)*velocity(i,j);
        }
        quantity_(i,0) *= mass(i,0);
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class FluctuatingVelocity
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  FluctuatingVelocity::FluctuatingVelocity(ATC_Method * atc,
    PerAtomQuantity<double> * atomVelocities,
    PerAtomQuantity<double> * atomMeanVelocities,
    AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,3,atomType),
    atomVelocities_(atomVelocities),
    atomMeanVelocities_(atomMeanVelocities)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomVelocities_) {
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(
        LammpsInterface::ATOM_VELOCITY, atomType);
    }
    if (!atomMeanVelocities_) {
      atomMeanVelocities_ = interscaleManager.per_atom_quantity(field_to_prolongation_name(VELOCITY));
    }

    atomVelocities_->register_dependence(this);
    atomMeanVelocities_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  FluctuatingVelocity::~FluctuatingVelocity()
  {
    atomVelocities_->remove_dependence(this);
    atomMeanVelocities_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void FluctuatingVelocity::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & velocity(atomVelocities_->quantity());
      const DENS_MAT & meanVelocity(atomMeanVelocities_->quantity());
      
      // q = m * (v dot v)
      for (int i = 0; i < quantity_.nRows(); i++) {
        for (int j = 0; j < velocity.nCols(); j++) {
          quantity_(i,j) = velocity(i,j) - meanVelocity(i,j);
        }
      }
    }
  }
  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class ChargeVelocity
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  ChargeVelocity::ChargeVelocity(ATC_Method * atc,
    PerAtomQuantity<double> * atomVelocities,
    FundamentalAtomQuantity * atomCharge,
    AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,3,atomType),
    fluctuatingVelocities_(atomVelocities),
    atomCharge_(atomCharge)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!fluctuatingVelocities_) {
      fluctuatingVelocities_ = interscaleManager.per_atom_quantity("AtomicFluctuatingVelocity");
    }
    if (!atomCharge_) {
      atomCharge_ = interscaleManager.fundamental_atom_quantity(
        LammpsInterface::ATOM_CHARGE, atomType);
    }
    fluctuatingVelocities_->register_dependence(this);
    atomCharge_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  ChargeVelocity::~ChargeVelocity()
  {
    fluctuatingVelocities_->remove_dependence(this);
    atomCharge_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void ChargeVelocity::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & velocity(fluctuatingVelocities_->quantity());
      const DENS_MAT & charge(atomCharge_->quantity());
      for (int i = 0; i < quantity_.nRows(); i++) {
        for (int j = 0; j < velocity.nCols(); j++) {
          quantity_(i,j) = charge(i,0)*velocity(i,j);
        }
      }
    }
  }
  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class SpeciesVelocity
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  SpeciesVelocity::SpeciesVelocity(ATC_Method * atc,
    PerAtomQuantity<double> * fluctuatingVelocities,
    PerAtomQuantity<double> * atomTypeVector,
    AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,3*(atomTypeVector->nCols()),atomType),
    fluctuatingVelocities_(fluctuatingVelocities),
    atomTypeVector_(atomTypeVector)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!fluctuatingVelocities_) {
      fluctuatingVelocities_ = interscaleManager.per_atom_quantity("AtomicFluctuatingVelocity");
    }
    if (!atomTypeVector_) {
      atomTypeVector_ = interscaleManager.per_atom_quantity("AtomicTypeVector");
    }
    fluctuatingVelocities_->register_dependence(this);
    atomTypeVector_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  SpeciesVelocity::~SpeciesVelocity()
  {
    fluctuatingVelocities_->remove_dependence(this);
    atomTypeVector_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void SpeciesVelocity::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & velocity(fluctuatingVelocities_->quantity());
      const DENS_MAT & tv(atomTypeVector_->quantity());
      
      for (int i = 0; i < quantity_.nRows(); i++) {
        int m = 0;
        for (int j = 0; j < velocity.nCols(); j++) {
          for (int k = 0; j < tv.nCols(); j++) {
            quantity_(i,m++) = tv(i,k)*velocity(i,j);
          }
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class TwiceFluctuatingKineticEnergy
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  TwiceFluctuatingKineticEnergy::TwiceFluctuatingKineticEnergy(ATC_Method * atc,
    PerAtomQuantity<double> * atomVelocities,
    PerAtomQuantity<double> * atomMasses,
    PerAtomQuantity<double> * atomMeanVelocities,
    AtomType atomType) :
    AtomicEnergyForTemperature(atc,atomType),
    atomVelocities_(atomVelocities),
    atomMasses_(atomMasses),
    atomMeanVelocities_(atomMeanVelocities)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomVelocities_) {
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY,
                                                                    atomType);
    }
    if (!atomMasses_) {
      atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                                atomType);
    }
    if (!atomMeanVelocities_) {
      atomMeanVelocities_ = interscaleManager.per_atom_quantity(field_to_prolongation_name(VELOCITY));
    }

    atomVelocities_->register_dependence(this);
    atomMasses_->register_dependence(this);
    atomMeanVelocities_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  TwiceFluctuatingKineticEnergy::~TwiceFluctuatingKineticEnergy()
  {
    atomVelocities_->remove_dependence(this);
    atomMasses_->remove_dependence(this);
    atomMeanVelocities_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void TwiceFluctuatingKineticEnergy::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & mass(atomMasses_->quantity());
      const DENS_MAT & velocity(atomVelocities_->quantity());
      const DENS_MAT & meanVelocity(atomMeanVelocities_->quantity());
      
      // q = m * (v dot v)
      double vRel;
      for (int i = 0; i < quantity_.nRows(); i++) {
        vRel = velocity(i,0) - meanVelocity(i,0);
        quantity_(i,0) = vRel*vRel;
        for (int j = 1; j < velocity.nCols(); j++) {
          vRel = velocity(i,j) - meanVelocity(i,j);
          quantity_(i,0) += vRel*vRel;;
        }
        quantity_(i,0) *= mass(i,0);
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class KineticTensor
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  KineticTensor::KineticTensor(ATC_Method * atc,
     PerAtomQuantity<double> * atomVelocities,
     PerAtomQuantity<double> * atomMasses,
     AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,6,atomType),
    atomVelocities_(atomVelocities),
    atomMasses_(atomMasses)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomVelocities_) {
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY,
                                                                    atomType);
    }
    if (!atomMasses_) {
      atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                                atomType);
    }

    atomVelocities_->register_dependence(this);
    atomMasses_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  KineticTensor::~KineticTensor()
  {
    atomVelocities_->remove_dependence(this);
    atomMasses_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void KineticTensor::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & mass(atomMasses_->quantity());
      const DENS_MAT & velocity(atomVelocities_->quantity());
      
      // K = m * (v \otimes v)
      for (int i = 0; i < quantity_.nRows(); i++) {
        double m = mass(i,0);
        double v[3] = {velocity(i,0),velocity(i,1),velocity(i,2)};
        quantity_(i,0) -= m*v[0]*v[0];
        quantity_(i,1) -= m*v[1]*v[1];
        quantity_(i,2) -= m*v[2]*v[2];
        quantity_(i,3) -= m*v[0]*v[1];
        quantity_(i,4) -= m*v[0]*v[2];
        quantity_(i,5) -= m*v[1]*v[2];
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class FluctuatingKineticTensor
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  FluctuatingKineticTensor::FluctuatingKineticTensor(ATC_Method * atc,
     PerAtomQuantity<double> * atomVelocities,
     PerAtomQuantity<double> * atomMasses,
     PerAtomQuantity<double> * atomMeanVelocities,
     AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,6,atomType),
    atomVelocities_(atomVelocities),
    atomMasses_(atomMasses)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomVelocities_) {
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY,
                                                                    atomType);
    }
    if (!atomMasses_) {
      atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                                atomType);
    }
    if (!atomMeanVelocities_) {
      atomMeanVelocities_ = interscaleManager.per_atom_quantity(field_to_prolongation_name(VELOCITY));
    }

    atomVelocities_->register_dependence(this);
    atomMasses_->register_dependence(this);
    atomMeanVelocities_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  FluctuatingKineticTensor::~FluctuatingKineticTensor()
  {
    atomVelocities_->remove_dependence(this);
    atomMasses_->remove_dependence(this);
    atomMeanVelocities_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void FluctuatingKineticTensor::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & mass(atomMasses_->quantity());
      const DENS_MAT & velocity(atomVelocities_->quantity());
      const DENS_MAT & meanVelocity(atomMeanVelocities_->quantity());
      
      // K = m * (v \otimes v)
      for (int i = 0; i < quantity_.nRows(); i++) {
        double m = mass(i,0);
        double v[3] = {velocity(i,0)-meanVelocity(i,0),
                       velocity(i,1)-meanVelocity(i,0),
                       velocity(i,2)-meanVelocity(i,0)};
        quantity_(i,0) -= m*v[0]*v[0];
        quantity_(i,1) -= m*v[1]*v[1];
        quantity_(i,2) -= m*v[2]*v[2];
        quantity_(i,3) -= m*v[0]*v[1];
        quantity_(i,4) -= m*v[0]*v[2];
        quantity_(i,5) -= m*v[1]*v[2];
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class MixedKePeEnergy
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  MixedKePeEnergy::MixedKePeEnergy(ATC_Method * atc,
                                   double keMultiplier,
                                   double peMultiplier,
                                   PerAtomQuantity<double> * twiceKineticEnergy,
                                   PerAtomQuantity<double> * potentialEnergy,
                                   AtomType atomType) :
    AtomicEnergyForTemperature(atc,atomType),
    keMultiplier_(keMultiplier),
    peMultiplier_(peMultiplier),
    twiceKineticEnergy_(twiceKineticEnergy),
    potentialEnergy_(potentialEnergy)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!twiceKineticEnergy_) {
       twiceKineticEnergy_ = interscaleManager.per_atom_quantity("AtomicTwiceKineticEnergy");
    }
    if (!potentialEnergy_) {
      potentialEnergy_ = interscaleManager.per_atom_quantity("AtomicFluctuatingPotentialEnergy");
    }

    twiceKineticEnergy_->register_dependence(this);
    potentialEnergy_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  MixedKePeEnergy::~MixedKePeEnergy()
  {
    twiceKineticEnergy_->remove_dependence(this);
    potentialEnergy_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void MixedKePeEnergy::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & twoKe(twiceKineticEnergy_->quantity());
      const DENS_MAT & pe(potentialEnergy_->quantity());

      // q = peScale * pe + keScale * ke
      quantity_ = pe;
      quantity_ *= peMultiplier_;
      quantity_ += (keMultiplier_/2.)*twoKe;
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class TotalEnergy
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  TotalEnergy::TotalEnergy(ATC_Method * atc,
                           PerAtomQuantity<double> * twiceKineticEnergy,
                           PerAtomQuantity<double> * potentialEnergy,
                           AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,atomType),
    twiceKineticEnergy_(twiceKineticEnergy),
    potentialEnergy_(potentialEnergy)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());

    if (!twiceKineticEnergy_) {
       twiceKineticEnergy_ = interscaleManager.per_atom_quantity("AtomicTwiceKineticEnergy");
    }
    if (!potentialEnergy_) {
      potentialEnergy_ = interscaleManager.per_atom_quantity("AtomicPotentialEnergy");
    }
    twiceKineticEnergy_->register_dependence(this);
    potentialEnergy_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  TotalEnergy::~TotalEnergy()
  {
    twiceKineticEnergy_->remove_dependence(this);
    potentialEnergy_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void TotalEnergy::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & twoKe(twiceKineticEnergy_->quantity());
      const DENS_MAT & pe(potentialEnergy_->quantity());
      quantity_ = pe;
      quantity_ += (0.5)*twoKe;
    }
  }
  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class FluctuatingPotentialEnergy
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  FluctuatingPotentialEnergy::FluctuatingPotentialEnergy(ATC_Method * atc,
                                                         PerAtomQuantity<double> * potentialEnergy,
                                                         PerAtomQuantity<double> * referencePotential,
                                                         AtomType atomType) :
    AtomicEnergyForTemperature(atc,atomType),
    potentialEnergy_(potentialEnergy),
    referencePotential_(referencePotential)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!potentialEnergy_) {
      potentialEnergy_ = interscaleManager.per_atom_quantity("AtomicPotentialEnergy");
    }
    if (!referencePotential_) {
       referencePotential_ = interscaleManager.per_atom_quantity("AtomicReferencePotential");
    }
    
    potentialEnergy_->register_dependence(this);
    referencePotential_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  FluctuatingPotentialEnergy::~FluctuatingPotentialEnergy()
  {
    potentialEnergy_->remove_dependence(this);
    referencePotential_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void FluctuatingPotentialEnergy::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & pe(potentialEnergy_->quantity());
      const DENS_MAT & refPe(referencePotential_->quantity());

      quantity_ = pe;
      quantity_ -= refPe;
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class DotTwiceKineticEnergy
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  DotTwiceKineticEnergy::DotTwiceKineticEnergy(ATC_Method * atc,
                                               PerAtomQuantity<double> * atomForces,
                                               PerAtomQuantity<double> * atomVelocities,
                                               AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,1,atomType),
    atomForces_(atomForces),
    atomVelocities_(atomVelocities)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomForces_)
      atomForces_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_FORCE,
                                                                    atomType);
    if (!atomVelocities_)
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY,
                                                                        atomType);

    atomForces_->register_dependence(this);
    atomVelocities_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  DotTwiceKineticEnergy::~DotTwiceKineticEnergy()
  {
    atomForces_->remove_dependence(this);
    atomVelocities_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void DotTwiceKineticEnergy::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & velocity(atomVelocities_->quantity());
      const DENS_MAT & force(atomForces_->quantity());

      for (int i = 0; i < quantity_.nRows(); i++) {
        quantity_(i,0) = velocity(i,0)*force(i,0);
        for (int j = 1; j < velocity.nCols(); j++) {
          quantity_(i,0) += velocity(i,j)*force(i,j);
        }
        quantity_(i,0) *= 2.;
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class VelocitySquared
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  VelocitySquared::VelocitySquared(ATC_Method * atc,
                                   PerAtomQuantity<double> * atomVelocities,
                                   AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,1,atomType),
    atomVelocities_(atomVelocities)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomVelocities_)
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY,
                                                                    atomType);
    atomVelocities_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  VelocitySquared::~VelocitySquared()
  {
    atomVelocities_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void VelocitySquared::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & velocity(atomVelocities_->quantity());

      for (int i = 0; i < quantity_.nRows(); i++) {
        quantity_(i,0) = velocity(i,0)*velocity(i,0);
        for (int j = 1; j < velocity.nCols(); j++) {
          quantity_(i,0) += velocity(i,j)*velocity(i,j);
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class LambdaSquared
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  LambdaSquared::LambdaSquared(ATC_Method * atc,
                               PerAtomQuantity<double> * atomMasses,
                               PerAtomQuantity<double> * atomVelocitiesSquared,
                               PerAtomQuantity<double> * atomLambdas,
                               AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,1,atomType),
    atomMasses_(atomMasses),
    atomVelocitiesSquared_(atomVelocitiesSquared),
    atomLambdas_(atomLambdas)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomMasses_) {
      atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                                    atomType);
    }
    atomMasses_->register_dependence(this);

    if (!atomVelocitiesSquared) {
      atomVelocitiesSquared_ = interscaleManager.per_atom_quantity("LambdaSquared");
    }
    atomVelocitiesSquared_->register_dependence(this);

    if (!atomLambdas_) {
      atomLambdas_ = interscaleManager.per_atom_quantity("AtomLambdaEnergy");
    }
    atomLambdas_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  LambdaSquared::~LambdaSquared()
  {
    atomMasses_->remove_dependence(this);
    atomVelocitiesSquared_->remove_dependence(this);
    atomLambdas_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void LambdaSquared::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & mass(atomMasses_->quantity());
      const DENS_MAT & velocitySquared(atomVelocitiesSquared_->quantity());
      const DENS_MAT & lambda(atomLambdas_->quantity());

      for (int i = 0; i < quantity_.nRows(); i++) {
        quantity_(i,0) = lambda(i,0)*lambda(i,0)*velocitySquared(i,0)/mass(i,0);
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomToType
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomToType::AtomToType(ATC_Method * atc,
                         int type,
                         AtomType atomType) :
    LargeToSmallAtomMap(atc,atomType),
    type_(type)
  {
    // DO NOTHING
    
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomToType::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<int>::reset();
      quantity_ = -1.;
      size_ = 0;
      const int * type = lammpsInterface_->atom_type();

      const Array<int> & quantityToLammps = atc_.atc_to_lammps_map();
      for (int i = 0; i < quantity_.nRows(); ++i) {
        int atomIdx = quantityToLammps(i);
        if (type[atomIdx] == type_) {
          quantity_(i,0) = size_;
          ++size_;
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomToGroup
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomToGroup::AtomToGroup(ATC_Method * atc,
                           int group,
                           AtomType atomType) :
    LargeToSmallAtomMap(atc,atomType),
    group_(group)
  {
    // DO NOTHING
    
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomToGroup::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<int>::reset();
      quantity_ = -1.;
      size_ = 0;
      const int * mask = lammpsInterface_->atom_mask();

      const Array<int> & quantityToLammps = atc_.atc_to_lammps_map();
      for (int i = 0; i < quantity_.nRows(); ++i) {
        int atomIdx = quantityToLammps(i);
        if (mask[atomIdx] & group_) {
          quantity_(i,0) = size_;
          ++size_;
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomToNodeset
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomToNodeset::AtomToNodeset(ATC_Method * atc,
                               SetDependencyManager<int> * subsetNodes,
                               PerAtomQuantity<int> * atomElement,
                               AtomType atomType) :
    LargeToSmallAtomMap(atc,atomType),
    subsetNodes_(subsetNodes),
    atomElement_(atomElement),
    feMesh_((atc->fe_engine())->fe_mesh())
  {
    if (!atomElement_) {
      atomElement_ = (atc->interscale_manager()).per_atom_int_quantity("AtomElement");
    }
    if (atomElement_) {
      atomElement_->register_dependence(this);
    }
    else {
      throw ATC_Error("PerAtomQuantityLibrary::AtomToRegulated - No AtomElement provided");
    }

    subsetNodes_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  quantity
  //--------------------------------------------------------
  void AtomToNodeset::reset() const
  {
    
    //so it has been commented out.
    /*
    if (need_reset()) {
      PerAtomQuantity<int>::reset();
      const INT_ARRAY & atomElement(atomElement_->quantity());
      const set<int> & subsetNodes(subsetNodes_->quantity());
      int nLocal = atomElement.nRows();
      quantity_.resize(nLocal,1);
      quantity_ = -1;
      size_ = 0;

      for (int i = 0; i < nLocal; ++i) {
        feMesh_->element_connectivity_unique(atomElement(i,0),_nodes_);
        for (int j = 0; j < _nodes_.size(); ++j) {
          if (subsetNodes.find(_nodes_(j)) != subsetNodes.end()) {
            quantity_(i,0) = size_;
            size_++;
            break;
          }
        }
      }
    }
    */
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomToElementset
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  AtomToElementset::AtomToElementset(ATC_Method * atc,
                                     MatrixDependencyManager<DenseMatrix, bool> * elementMask,
                                     PerAtomQuantity<int> * atomElement,
                                     AtomType atomType) :
    LargeToSmallAtomMap(atc,atomType),
    elementMask_(elementMask),
    atomElement_(atomElement),
    feMesh_((atc->fe_engine())->fe_mesh())
  {
    if (!atomElement_) {
      atomElement_ = (atc->interscale_manager()).per_atom_int_quantity("AtomElement");
    }
    if (atomElement_) {
      atomElement_->register_dependence(this);
    }
    else {
      throw ATC_Error("PerAtomQuantityLibrary::AtomToRegulated - No AtomElement provided");
    }

    elementMask_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomToElementset::~AtomToElementset()
  {
    atomElement_->remove_dependence(this);
    elementMask_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  quantity
  //--------------------------------------------------------
  void AtomToElementset::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<int>::reset();
      const INT_ARRAY & atomElement(atomElement_->quantity());
      const DenseMatrix<bool> & elementMask(elementMask_->quantity());
      int nLocal = atomElement.nRows();
      quantity_.resize(nLocal,1);
      quantity_ = -1;
      size_ = 0;

      for (int i = 0; i < nLocal; ++i) {
        if (elementMask(atomElement(i,0),0)) {
          quantity_(i,0) = size_;
          size_++;
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class MappedAtomQuantity
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  MappedAtomQuantity::MappedAtomQuantity(ATC_Method * atc,
                                         PerAtomQuantity<double> * source,
                                         LargeToSmallAtomMap * map,
                                         AtomType atomType) :
    ProtectedMappedAtomQuantity<double>(atc,map,source->nCols(),atomType),
    source_(source),
    map_(map)
  {
    source_->register_dependence(this);
    map_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  reset_quantity
  //--------------------------------------------------------
  void MappedAtomQuantity::reset() const
  {
    if (needReset_) {
      ProtectedMappedAtomQuantity<double>::reset();
      const DENS_MAT & source(source_->quantity());
      const INT_ARRAY & map(map_->quantity());
      quantity_.resize(map_->size(),source.nCols());
      
      for (int i = 0; i < source.nRows(); i++) {
        int idx = map(i,0);
        if (idx > -1) {
          for (int j = 0; j < source.nCols(); j++) {
            quantity_(idx,j) = source(i,j);
          }
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class VelocitySquaredMapped
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  VelocitySquaredMapped::VelocitySquaredMapped(ATC_Method * atc,
                                               MatrixDependencyManager<DenseMatrix, int> * atomMap,
                                               PerAtomQuantity<double> * atomVelocities,
                                               AtomType atomType) :
    ProtectedMappedAtomQuantity<double>(atc,atomMap,1,atomType),
    atomVelocities_(atomVelocities)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomVelocities_)
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY,
                                                                        atomType);
    atomVelocities_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  VelocitySquaredMapped::~VelocitySquaredMapped()
  {
    atomVelocities_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void VelocitySquaredMapped::reset() const
  {
    if (need_reset()) {
      ProtectedMappedAtomQuantity<double>::reset();
      const DENS_MAT & velocity(atomVelocities_->quantity());
      const INT_ARRAY & atomMap(atomMap_->quantity());

      for (int i = 0; i < atomMap.nRows(); i++) {
        int idx = atomMap(i,0);
        if (idx > -1) {
          quantity_(idx,0) = velocity(i,0)*velocity(i,0);
          for (int j = 1; j < velocity.nCols(); j++) {
            quantity_(idx,0) += velocity(i,j)*velocity(i,j);
          }
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class LambdaSquaredMapped
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  LambdaSquaredMapped::LambdaSquaredMapped(ATC_Method * atc,
                                           MatrixDependencyManager<DenseMatrix, int> * atomMap,
                                           PerAtomQuantity<double> * atomMasses,
                                           PerAtomQuantity<double> * atomVelocitiesSquared,
                                           PerAtomQuantity<double> * atomLambdas,
                                           AtomType atomType) :
    ProtectedMappedAtomQuantity<double>(atc,atomMap,1,atomType),
    atomMasses_(atomMasses),
    atomVelocitiesSquared_(atomVelocitiesSquared),
    atomLambdas_(atomLambdas)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomMasses_) {
      atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                                    atomType);
    }
    atomMasses_->register_dependence(this);

    if (!atomVelocitiesSquared) {
      atomVelocitiesSquared_ = interscaleManager.per_atom_quantity("LambdaSquared");
    }
    atomVelocitiesSquared_->register_dependence(this);

    if (!atomLambdas_) {
      atomLambdas_ = interscaleManager.per_atom_quantity("AtomLambdaEnergy");
    }
    atomLambdas_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  LambdaSquaredMapped::~LambdaSquaredMapped()
  {
    atomMasses_->remove_dependence(this);
    atomVelocitiesSquared_->remove_dependence(this);
    atomLambdas_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void LambdaSquaredMapped::reset() const
  {
    if (need_reset()) {
      ProtectedMappedAtomQuantity<double>::reset();
      const DENS_MAT & mass(atomMasses_->quantity());
      const DENS_MAT & velocitySquared(atomVelocitiesSquared_->quantity());
      const DENS_MAT & lambda(atomLambdas_->quantity());
      const INT_ARRAY & atomMap(atomMap_->quantity());

      for (int i = 0; i < atomMap.nRows(); i++) {
        int idx = atomMap(i,0);
        if (idx > -1) {
          quantity_(idx,0) = lambda(i,0)*lambda(i,0)*velocitySquared(idx,0)/mass(i,0);
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class HeatCapacity
  //    computes the classical atomic heat capacity
  //--------------------------------------------------------
  //--------------------------------------------------------
  
  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  HeatCapacity::HeatCapacity(ATC_Method * atc, AtomType atomType) :
    ConstantQuantity<double>(atc,1,1,atomType)
  {
    constant_ = (atc->nsd())*(lammpsInterface_->kBoltzmann());
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomicVelocityRescaleFactor
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomicVelocityRescaleFactor::AtomicVelocityRescaleFactor(ATC_Method * atc,
                                                           PerAtomQuantity<double> * atomLambdas,
                                                           AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,1,atomType),
    atomLambdas_(atomLambdas)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomLambdas) {
      atomLambdas_ = interscaleManager.per_atom_quantity("AtomLambdaEnergy");
    }
    atomLambdas_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomicVelocityRescaleFactor::~AtomicVelocityRescaleFactor()
  {
    atomLambdas_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomicVelocityRescaleFactor::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & lambda(atomLambdas_->quantity());

      for (int i = 0; i < quantity_.nRows(); i++) {
        quantity_(i,0) = sqrt(lambda(i,0));
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomicFluctuatingVelocityRescaled
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomicFluctuatingVelocityRescaled::AtomicFluctuatingVelocityRescaled(ATC_Method * atc,
                                                                       PerAtomQuantity<double> * atomRescaleFactor,
                                                                       PerAtomQuantity<double> * atomFluctuatingVelocity,
                                                                       AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,atc->nsd(),atomType),
    atomRescaleFactor_(atomRescaleFactor),
    atomFluctuatingVelocity_(atomFluctuatingVelocity)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomRescaleFactor_) {
      atomRescaleFactor_ = interscaleManager.per_atom_quantity("AtomVelocityRescaling");
    }
    if (!atomFluctuatingVelocity_) {
      atomFluctuatingVelocity_ = interscaleManager.per_atom_quantity("AtomicFluctuatingVelocity");
    }

    atomRescaleFactor_->register_dependence(this);
    atomFluctuatingVelocity_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomicFluctuatingVelocityRescaled::~AtomicFluctuatingVelocityRescaled()
  {
    atomRescaleFactor_->remove_dependence(this);
    atomFluctuatingVelocity_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomicFluctuatingVelocityRescaled::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & factor(atomRescaleFactor_->quantity());
      const DENS_MAT & v(atomFluctuatingVelocity_->quantity());

      for (int i = 0; i < quantity_.nRows(); i++) {
        for (int j = 0; j < quantity_.nCols(); j++) {
          quantity_(i,j) = factor(i,0)*v(i,j);
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomicCombinedRescaleThermostatError
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomicCombinedRescaleThermostatError::AtomicCombinedRescaleThermostatError(ATC_Method * atc,
                                                                             PerAtomQuantity<double> * atomFluctuatingMomentumRescaled,
                                                                             PerAtomQuantity<double> * atomMeanVelocity,
                                                                             PerAtomQuantity<double> * atomStreamingVelocity,
                                                                             PerAtomQuantity<double> * atomMass,
                                                                             AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,1,atomType),
    atomFluctuatingMomentumRescaled_(atomFluctuatingMomentumRescaled),
    atomMeanVelocity_(atomMeanVelocity),
    atomStreamingVelocity_(atomStreamingVelocity),
    atomMass_(atomMass)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomFluctuatingMomentumRescaled_) {
      atomFluctuatingMomentumRescaled_ = interscaleManager.per_atom_quantity("AtomFluctuatingMomentumRescaled");
    }
    if (!atomMeanVelocity_) {
      atomMeanVelocity_ = interscaleManager.per_atom_quantity(field_to_prolongation_name(VELOCITY));
    }
    if (!atomStreamingVelocity_) {
      atomStreamingVelocity_ = interscaleManager.per_atom_quantity("AtomLambdaMomentum");
    }
    if (!atomMass_) {
      atomMass_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                              atomType);
    }

    atomFluctuatingMomentumRescaled_->register_dependence(this);
    atomMeanVelocity_->register_dependence(this);
    atomStreamingVelocity_->register_dependence(this);
    atomMass_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomicCombinedRescaleThermostatError::~AtomicCombinedRescaleThermostatError()
  {
    atomFluctuatingMomentumRescaled_->remove_dependence(this);
    atomMeanVelocity_->remove_dependence(this);
    atomStreamingVelocity_->remove_dependence(this);
    atomMass_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomicCombinedRescaleThermostatError::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & m(atomMass_->quantity());
      const DENS_MAT & p(atomFluctuatingMomentumRescaled_->quantity());
      const DENS_MAT & v(atomMeanVelocity_->quantity());
      const DENS_MAT & s(atomStreamingVelocity_->quantity());

      double diff;
      for (int i = 0; i < quantity_.nRows(); i++) {
        diff = s(i,0)-v(i,0);
        quantity_(i,0) = 2.*p(i,0)*diff + m(i,0)*diff*diff;
        for (int j = 1; j < p.nCols(); j++) {
          diff = s(i,j)-v(i,j);
          quantity_(i,0) += 2.*p(i,j)*diff + m(i,0)*diff*diff;
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomicThermostatForce
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomicThermostatForce::AtomicThermostatForce(ATC_Method * atc,
                                               PerAtomQuantity<double> * atomLambdas,
                                               PerAtomQuantity<double> * atomVelocities,
                                               AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,atc->nsd(),atomType),
    atomLambdas_(atomLambdas),
    atomVelocities_(atomVelocities)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomLambdas) {
      atomLambdas_ = interscaleManager.per_atom_quantity("AtomLambdaEnergy");
    }
    if (!atomVelocities_) {
      atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY,
                                                                    atomType);
    }
    
    atomLambdas_->register_dependence(this);
    atomVelocities_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomicThermostatForce::~AtomicThermostatForce()
  {
    atomLambdas_->remove_dependence(this);
    atomVelocities_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomicThermostatForce::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & v(atomVelocities_->quantity());
      const DENS_MAT & lambda(atomLambdas_->quantity());
      
      // force = -lambda*v
      quantity_ = v;
      quantity_ *= lambda;
      quantity_ *= -1.;
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomicKinetostatForceDisplacement
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomicKinetostatForceDisplacement::AtomicKinetostatForceDisplacement(ATC_Method * atc,
                                               PerAtomQuantity<double> * atomLambda,
                                               PerAtomQuantity<double> * atomMass,
                                               AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,atc->nsd(),atomType),
    atomLambda_(atomLambda),
    atomMass_(atomMass)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomLambda) {
      atomLambda_ = interscaleManager.per_atom_quantity("AtomLambdaMomentum");
    }
    if (!atomMass_) {
      atomMass_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS,
                                                              atomType);
    }
    
    atomLambda_->register_dependence(this);
    atomMass_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomicKinetostatForceDisplacement::~AtomicKinetostatForceDisplacement()
  {
    atomLambda_->remove_dependence(this);
    atomMass_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomicKinetostatForceDisplacement::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & m(atomMass_->quantity());
      const DENS_MAT & lambda(atomLambda_->quantity());
      double timeFactor = time_step_factor(0.5*atc_.dt());
      
      //force = -lambda*m*(timestep factor)
      quantity_ = lambda;
      quantity_ *= m;
      quantity_ *= -1.*timeFactor;
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class AtomicKinetostatForceStress
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  AtomicKinetostatForceStress::AtomicKinetostatForceStress(ATC_Method * atc,
                                                           PerAtomQuantity<double> * atomLambda,
                                                           AtomType atomType) :
    ProtectedAtomQuantity<double>(atc,atc->nsd(),atomType),
    atomLambda_(atomLambda)
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomLambda_) {
      atomLambda_ = interscaleManager.per_atom_quantity("AtomLambdaMomentum");
    }
    
    atomLambda_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  AtomicKinetostatForceStress::~AtomicKinetostatForceStress()
  {
    atomLambda_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void AtomicKinetostatForceStress::reset() const
  {
    if (need_reset()) {
      PerAtomQuantity<double>::reset();
      const DENS_MAT & lambda(atomLambda_->quantity());
      
      // force = -lambda
      quantity_ = lambda;
      quantity_ *= -1.;
    }
  }
  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class PerAtomKernelFunction
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  PerAtomKernelFunction::PerAtomKernelFunction(ATC_Method * atc,
                                               PerAtomQuantity<double> * atomPositions,
                                               AtomType atomType) :
    ProtectedAtomSparseMatrix<double>(atc,(atc->fe_engine())->num_nodes(),atc->accumulant_bandwidth(),atomType),
    atomPositions_(atomPositions),
    feEngine_(atc->fe_engine())
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomPositions_) {
      atomPositions_ = interscaleManager.per_atom_quantity("AtomicCoarseGrainingPositions");
    }
    atomPositions_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  PerAtomKernelFunction::~PerAtomKernelFunction()
  {
    atomPositions_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void PerAtomKernelFunction::reset() const
  {
    if (need_reset()) {
      PerAtomSparseMatrix<double>::reset();
      const DENS_MAT & positions(atomPositions_->quantity());
      if (positions.nRows() > 0) {
        feEngine_->evaluate_kernel_functions(positions,quantity_);
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class PerAtomShapeFunction
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  constructor
  //--------------------------------------------------------
  PerAtomShapeFunction::PerAtomShapeFunction(ATC_Method * atc,
                                             PerAtomQuantity<double> * atomPositions,
                                             PerAtomQuantity<int> * atomElements,
                                             AtomType atomType) :
    ProtectedAtomSparseMatrix<double>(atc,atc->num_nodes(),(atc->fe_engine())->num_nodes_per_element(),atomType),
    atomPositions_(atomPositions),
    atomElements_(atomElements),
    feEngine_(atc->fe_engine())
  {
    InterscaleManager & interscaleManager(atc->interscale_manager());
    if (!atomPositions_) {
      atomPositions_ = interscaleManager.per_atom_quantity("AtomicCoarseGrainingPositions");
    }
    if (!atomElements_) {
      atomElements_ = interscaleManager.per_atom_int_quantity("AtomElement");
    }

    atomPositions_->register_dependence(this);
    atomElements_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  destructor
  //--------------------------------------------------------
  PerAtomShapeFunction::~PerAtomShapeFunction()
  {
    atomPositions_->remove_dependence(this);
    atomElements_->remove_dependence(this);
  }

  //--------------------------------------------------------
  //  reset
  //--------------------------------------------------------
  void PerAtomShapeFunction::reset() const
  {
    if (need_reset()) {
      PerAtomSparseMatrix<double>::reset();
      const DENS_MAT & positions(atomPositions_->quantity());
      const INT_ARRAY & elements(atomElements_->quantity());
      if (positions.nRows() > 0) {
        feEngine_->evaluate_shape_functions(positions,
                                            elements,
                                            quantity_);
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class LambdaCouplingMatrix
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  LambdaCouplingMatrix::LambdaCouplingMatrix(ATC_Method * atc,
                                             MatrixDependencyManager<DenseMatrix, int> * nodeToOverlapMap,
                                             SPAR_MAN * shapeFunction) :
    ProtectedMappedAtomSparseMatrix<double>(atc,nodeToOverlapMap->size(),
                                            (atc->fe_engine())->num_nodes_per_element(),
                                            INTERNAL),
    shapeFunction_(shapeFunction),
    nodeToOverlapMap_(nodeToOverlapMap)
  {
    if (!shapeFunction_) {
      shapeFunction_ = (atc->interscale_manager()).per_atom_sparse_matrix("Interpolant");;
    }
    if (!nodeToOverlapMap_) {
      nodeToOverlapMap_ = (atc->interscale_manager()).dense_matrix_int("NodeToOverlapMap");
    }
    
    shapeFunction_->register_dependence(this);
    nodeToOverlapMap_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  reset_quantity
  //--------------------------------------------------------
  void LambdaCouplingMatrix::reset() const
  {
    if (need_reset()) {
      PerAtomSparseMatrix<double>::reset();
      int nNodeOverlap = nodeToOverlapMap_->size();
      const SPAR_MAT & shapeFunction(shapeFunction_->quantity());
      quantity_.reset(shapeFunction.nRows(),nNodeOverlap); // number of atoms X number of nodes overlapping MD region
      const INT_ARRAY nodeToOverlapMap(nodeToOverlapMap_->quantity());
    
      for (int i = 0; i < shapeFunction.size(); ++i) {
        TRIPLET<double> myTriplet = shapeFunction.triplet(i);
        int myCol = nodeToOverlapMap(myTriplet.j,0);
        if (myCol > -1) {
          quantity_.set(myTriplet.i,myCol,myTriplet.v);
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class LocalLambdaCouplingMatrix
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  LocalLambdaCouplingMatrix::LocalLambdaCouplingMatrix(ATC_Method * atc,
                                                       MatrixDependencyManager<DenseMatrix, int> * lambdaAtomMap,
                                                       MatrixDependencyManager<DenseMatrix, int> * nodeToOverlapMap,
                                                       SPAR_MAN * shapeFunction) :
    LambdaCouplingMatrix(atc,nodeToOverlapMap,shapeFunction),
    lambdaAtomMap_(lambdaAtomMap)
  {
    if (!lambdaAtomMap_) {
      lambdaAtomMap_ = (atc->interscale_manager()).dense_matrix_int("LambdaAtomMap");
    }
    
    lambdaAtomMap_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  reset_quantity
  //--------------------------------------------------------
  void LocalLambdaCouplingMatrix::reset() const
  {
    if (need_reset()) {
      PerAtomSparseMatrix<double>::reset();
      int nNodeOverlap = nodeToOverlapMap_->size();
      int nLocalLambda = lambdaAtomMap_->size();

      quantity_.reset(nLocalLambda,nNodeOverlap); // number of regulated atoms X number of nodes containing them
      const SPAR_MAT & shapeFunction(shapeFunction_->quantity());
      const INT_ARRAY nodeToOverlapMap(nodeToOverlapMap_->quantity());
      const INT_ARRAY lambdaAtomMap(lambdaAtomMap_->quantity());
      
      for (int i = 0; i < shapeFunction.size(); ++i) {
        TRIPLET<double> myTriplet = shapeFunction.triplet(i);
        int myRow = lambdaAtomMap(myTriplet.i,0);
        int myCol = nodeToOverlapMap(myTriplet.j,0);
        if ((myRow > -1) && (myCol > -1)) {
          quantity_.set(myRow,myCol,myTriplet.v);
        }
      }
    }
  }

  //--------------------------------------------------------
  //--------------------------------------------------------
  //  Class GhostCouplingMatrix
  //--------------------------------------------------------
  //--------------------------------------------------------

  //--------------------------------------------------------
  //  Constructor
  //--------------------------------------------------------
  GhostCouplingMatrix::GhostCouplingMatrix(ATC_Method * atc,
                                           SPAR_MAN * shapeFunction,
                                           SetDependencyManager<int> * subsetNodes,
                                           MatrixDependencyManager<DenseMatrix, int> * nodeToOverlapMap) :
    LambdaCouplingMatrix(atc,nodeToOverlapMap,shapeFunction),
    subsetNodes_(subsetNodes)
  {
    subsetNodes_->register_dependence(this);
  }

  //--------------------------------------------------------
  //  reset_quantity
  //--------------------------------------------------------
  void GhostCouplingMatrix::reset() const
  {
    if (need_reset()) {
      PerAtomSparseMatrix<double>::reset();
      const SPAR_MAT & shapeFunction(shapeFunction_->quantity());
      const set<int> & subsetNodes(subsetNodes_->quantity());
      int nNodes = nodeToOverlapMap_->nRows();
      int nLocalGhost = shapeFunction.nRows();
      quantity_.reset(nLocalGhost,nNodes);
      const INT_ARRAY nodeToOverlapMap(nodeToOverlapMap_->quantity());
      for (int i = 0; i < shapeFunction.size(); ++i) {
        TRIPLET<double> myTriplet = shapeFunction.triplet(i);
        int myCol = myTriplet.j;
        if (nodeToOverlapMap(myCol,0) > -1) {
          quantity_.set(myTriplet.i,myCol,myTriplet.v);
        }
      }
      quantity_.compress();
      
      //int nNodes = nodeToOverlapMap_->nRows();
      _activeNodes_.reset(nNodes);
      for (int i = 0; i < nNodes; ++i) {
        if (subsetNodes.find(i) != subsetNodes.end()) {
          _activeNodes_(i) = 1.;
        }
      }

      //const SPAR_MAT & shapeFunction(shapeFunction_->quantity());
      //int nLocalGhost = shapeFunction.nRows();
      _weights_ = shapeFunction*_activeNodes_;
      _weightMatrix_.resize(nLocalGhost,nLocalGhost);
      for (int i = 0; i < nLocalGhost; ++i) {
        _weightMatrix_(i,i) = 1./_weights_(i);
      }
      quantity_ = _weightMatrix_*quantity_;
    }
  }
};