diff --git a/src/md/lammps/domain_lammps.cc b/src/md/lammps/domain_lammps.cc
index 0056384..b91e3c7 100644
--- a/src/md/lammps/domain_lammps.cc
+++ b/src/md/lammps/domain_lammps.cc
@@ -1,1036 +1,1037 @@
 /**
  * @file   domain_lammps.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Jaehyun Cho <jaehyun.cho@epfl.ch>
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   Thu Jul 31 22:41:23 2014
  *
  * @brief  This is the model wrapping LAMMPS
  *
  * @section LICENSE
  *
  * Copyright INRIA and CEA
  *
  * The LibMultiScale is a C++ parallel framework for the multiscale
  * coupling methods dedicated to material simulations. This framework
  * provides an API which makes it possible to program coupled simulations
  * and integration of already existing codes.
  *
  * This Project was initiated in a collaboration between INRIA Futurs Bordeaux
  * within ScAlApplix team and CEA/DPTA Ile de France.
  * The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
  * within the LSMS/ENAC laboratory.
  *
  * This software is governed by the CeCILL-C license under French law and
  * abiding by the rules of distribution of free software.  You can  use,
  * modify and/ or redistribute the software under the terms of the CeCILL-C
  * license as circulated by CEA, CNRS and INRIA at the following URL
  * "http://www.cecill.info".
  *
  * As a counterpart to the access to the source code and  rights to copy,
  * modify and redistribute granted by the license, users are provided only
  * with a limited warranty  and the software's author,  the holder of the
  * economic rights,  and the successive licensors  have only  limited
  * liability.
  *
  * In this respect, the user's attention is drawn to the risks associated
  * with loading,  using,  modifying and/or developing or reproducing the
  * software by the user in light of its specific status of free software,
  * that may mean  that it is complicated to manipulate,  and  that  also
  * therefore means  that it is reserved for developers  and  experienced
  * professionals having in-depth computer knowledge. Users are therefore
  * encouraged to load and test the software's suitability as regards their
  * requirements in conditions enabling the security of their systems and/or
  * data to be ensured and,  more generally, to use and operate it in the
  * same conditions as regards security.
  *
  * The fact that you are presently reading this means that you have had
  * knowledge of the CeCILL-C license and that you accept its terms.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 //#define CHECK_STABILITY
 /* -------------------------------------------------------------------------- */
 #include "domain_lammps.hh"
 #include "factory_multiscale.hh"
 #include "import_lammps.hh"
 #include "lammps_common.hh"
 #include "lm_common.hh"
 #include "reference_manager_lammps.hh"
 #include "trace_atom.hh"
 /* -------------------------------------------------------------------------- */
 // LAMMPS includes
 #include <angle.h>
 #include <bond.h>
 #include <change_box.h>
 #include <dihedral.h>
 #include <fix_minimize.h>
 #include <improper.h>
 #include <integrate.h>
 #include <kspace.h>
 #include <min_cg.h>
 #include <pair.h>
 #include <variable.h>
 #include <verlet.h>
 //#include <timer.h>
 /* -------------------------------------------------------------------------- */
 #include <functional>
 #include <iomanip>
 #include <iostream>
 #include <iterator>
 /* -------------------------------------------------------------------------- */
 #include "compute_lammps.hh"
 #include "geometry_manager.hh"
 #include "lm_communicator.hh"
 #include "lm_parser.hh"
 /* -------------------------------------------------------------------------- */
 LAMMPS_NS::LAMMPS *lammps_main_object = nullptr;
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_LIBMULTISCALE__
 
 template <UInt Dim> void DomainLammps<Dim>::init() {
   DomainAtomic<ContainerLammps<Dim>, Dim>::init();
 
   if (!this->getCommGroup().amIinGroup())
     return;
 
   this->initModel();
   this->update->whichflag = 1;
   LAMMPS_NS::LAMMPS::init();
   this->update->integrate->setup(0);
   this->initDOFs();
   this->ref_manager->setState(true);
   this->atoms.getBoundingBox() = this->getDomainBoundingBox();
+  this->updateAcceleration();
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> ArrayView DomainLammps<Dim>::gradient() {
   return getField(_force);
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> ArrayView DomainLammps<Dim>::primal() {
   return getField(_position);
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim>
 auto DomainLammps<Dim>::primalTimeDerivative() -> ArrayView {
   return getField(_velocity);
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> auto DomainLammps<Dim>::mass() -> ArrayView {
   return getField(_mass);
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> auto DomainLammps<Dim>::acceleration() -> ArrayView {
   return getField(_acceleration);
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim>
 auto DomainLammps<Dim>::getField(FieldType field_type) -> ArrayView {
   switch (field_type) {
   case _position0:
     return ArrayView(&lammps_main_object->atom->x0[0][0],
                      lammps_main_object->atom->nlocal, 3);
   case _velocity:
     return ArrayView(&lammps_main_object->atom->v[0][0],
                      lammps_main_object->atom->nlocal, 3);
   case _position:
     return ArrayView(&lammps_main_object->atom->x[0][0],
                      lammps_main_object->atom->nlocal, 3);
   case _force:
     return ArrayView(&lammps_main_object->atom->f[0][0],
                      lammps_main_object->atom->nlocal, 3);
   case _mass:
     if (not is_mass_per_particle)
       return ArrayView(&lammps_main_object->atom->mass[0],
                        lammps_main_object->atom->ntypes, 1);
     return ArrayView(&lammps_main_object->atom->rmass[0],
                      lammps_main_object->atom->natoms, 1);
 
   case _acceleration:
     return ArrayView(this->acceleration_field.data(),
                      lammps_main_object->atom->nlocal, 3);
 
   case _angular_velocity:
     if (not lammps_main_object->atom->omega) {
       LM_FATAL("created atoms do not have angle degrees of freedom");
     }
     return ArrayView(&lammps_main_object->atom->omega[0][0],
                      lammps_main_object->atom->nlocal, 3);
 
   case _angular_acceleration:
     return ArrayView(this->angular_acceleration_field.data(),
                      lammps_main_object->atom->nlocal, 3);
   case _torque:
     if (not lammps_main_object->atom->torque) {
       LM_FATAL("created atoms do not have angle degrees of freedom");
     }
     return ArrayView(&lammps_main_object->atom->torque[0][0],
                      lammps_main_object->atom->nlocal, 3);
 
   default:
     LM_FATAL("Not accessible field: " << field_type);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim> void DomainLammps<Dim>::loadLammpsConfigFile() {
 
   // forward LM variables to lammps
 
   auto &vars = Parser::getAlgebraicVariables();
   for (auto &&[varname, val] : vars) {
     std::stringstream varval;
     varval << std::setprecision(16) << std::scientific << val;
     LAMMPS_NS::Variable *lmpvar = LAMMPS_NS::LAMMPS::input->variable;
 
     char cvarname[256];
     varname.copy(cvarname, 255);
     cvarname[varname.length()] = 0;
 
     char cvarval[30];
     varval.str().copy(cvarval, 29);
     cvarval[varval.str().length()] = 0;
     if (lmpvar->find(cvarname) >= 0)
       LM_FATAL("cannot export LibMultiscale variable "
                << varname << " to LAMMPS since variable already exists");
 
     DUMP(cvarname << " = " << cvarval, DBG_INFO);
     lmpvar->add_equal(cvarname, cvarval);
   }
 
   std::string file_to_read = lammpsfile;
   DUMP("lammpsfile"
            << " = " << lammpsfile,
        DBG_INFO);
 
   if (create_header_flag && this->getCommGroup().getMyRank() == 0) {
     std::string unit_name;
     if (code_unit_system == metal_unit_system)
       unit_name = "metal";
     else if (code_unit_system == real_unit_system)
       unit_name = "real";
     else if (code_unit_system == si_unit_system)
       unit_name = "si";
     else
       LM_FATAL("unknown units");
 
     file_to_read = "./generated_lammps.config";
     std::ofstream file_lammps(file_to_read.c_str());
 
     file_lammps << "units " << unit_name << std::endl;
     file_lammps << "dimension " << Dim << std::endl;
     file_lammps << "atom_style atomic" << std::endl;
     file_lammps << "boundary ";
     for (UInt i = 0; i < 3; ++i) {
       file_lammps << boundaries[i] << " ";
     }
     file_lammps << std::endl;
     file_lammps << "lattice " << lattice << " " << std::setprecision(15)
                 << static_cast<Real>(lattice_size);
 
     file_lammps << " orient x ";
     for (UInt i = 0; i < 3; ++i)
       file_lammps << lattice_orient_x[i] << " ";
     if (Dim > 1) {
       file_lammps << "orient y ";
       for (UInt i = 0; i < 3; ++i)
         file_lammps << lattice_orient_y[i] << " ";
     }
     if (Dim == 3) {
       file_lammps << "orient z ";
       for (UInt i = 0; i < 3; ++i)
         file_lammps << lattice_orient_z[i] << " ";
     }
 
     file_lammps << "spacing ";
     for (UInt i = 0; i < 3; ++i)
       file_lammps << lattice_spacing[i] << " ";
 
     file_lammps << "origin ";
     for (UInt i = 0; i < Dim; ++i)
       file_lammps << Real(lattice_origin[i]) << " ";
     for (UInt i = Dim; i < 3; ++i)
       file_lammps << 0.0 << " ";
     file_lammps << std::endl;
 
     file_lammps << "region box block ";
     for (UInt i = 0; i < 6; ++i)
       file_lammps << replica[i] << " ";
     file_lammps << std::endl;
     file_lammps << "create_box      1 box" << std::endl;
     file_lammps << "create_atoms    1 box" << std::endl;
     file_lammps << std::endl << std::endl;
 
     if (triclinic == 1 && floorf(tilt[1] * 100 + 0.5) / 100 != 0.0) {
       //      double value;
       std::string direction;
       Real tilting = floorf(tilt[1] * 100 + 0.5) / 100;
       if (tilt[0] == 0.0) {
         // value = tilt[1]/lattice_size/lattice_spacing[0];
         direction = "xy";
       } else if (tilt[0] == 1.0) {
         // value = tilt[1]/lattice_size/lattice_spacing[0];
         direction = "xz";
       } else if (tilt[0] == 2.0) {
         // value = tilt[1]/lattice_size/lattice_spacing[0];
         direction = "yz";
       } else
         LM_FATAL("tilt must be happend in xy(0), xz(1) or yz(2) axis!!!");
 
       file_lammps << "change_box      triclinic" << std::endl;
       file_lammps << "displace_box    all " << direction << " final " << tilting
                   << " remap none units box" << std::endl;
       // file_lammps << "displace_box    all "<<direction<<" final "<<value<<"
       // remap none"<< std::endl;
     }
 
     std::ifstream in(lammpsfile.c_str());
     std::istreambuf_iterator<char> src(in.rdbuf());
     std::istreambuf_iterator<char> end;
     std::ostream_iterator<char> dest(file_lammps);
     std::copy(src, end, dest);
   }
   DUMP("opening lammps config file " << file_to_read, DBG_INFO_STARTUP);
   infile = fopen(file_to_read.c_str(), "rb");
   if (infile == NULL)
     LM_FATAL("cannot open " << file_to_read << " as a lammps config file !!!");
   input->file();
   DUMP("lammps file loaded", DBG_INFO);
 
   // register lammps computes to libmultiscale system
   UInt ncompute = this->modify->ncompute;
 
   for (UInt i = 0; i < ncompute; ++i) {
     auto _ptr = std::make_shared<ComputeLammps<LAMMPS_NS::Compute>>(
         *(this->modify->compute[i]), *this);
     _ptr->acquireContext(*this);
     FilterManager::getManager().addObject(_ptr);
   }
 
   // register lammps computes to libmultiscale system
   UInt nfix = this->modify->nfix;
 
   for (UInt i = 0; i < nfix; ++i) {
     auto _ptr = std::make_shared<ComputeLammps<LAMMPS_NS::Fix>>(
         *(this->modify->fix[i]), *this);
     _ptr->acquireContext(*this);
     FilterManager::getManager().addObject(_ptr);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> void DomainLammps<Dim>::changeBox() {
   auto &geom_manager = GeometryManager::getManager();
   auto c = geom_manager.getGeometry(this->newGeomBox)->getBoundingBox();
   this->setDomainBoundingBox(c);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim> void DomainLammps<Dim>::initModel() {
 
   lammps_main_object = this;
 
   ref_manager = std::make_shared<ReferenceManagerLammps<Dim>>(this->atoms);
   this->atoms.setRefManager(ref_manager);
   ref_manager->acquireContext(*this);
 
   ImportLammpsInterface::createImportLammps(ref_manager);
 
   if (this->getGeom().getDim() == 1)
     LM_FATAL("LAMMPS cannot work in 1D");
 
   ImportLammpsInterface::setGeom(&this->getGeom());
 
   // for (auto &&_type : geom_by_type) {
   //   Geometry &g = GeometryManager::getManager().getObject(_type.second);
   //   UInt itype = _type.first;
   //   ImportLammpsInterface::setGeom(&g, itype);
   // }
 
   loadLammpsConfigFile();
 
   this->atoms.init();
 
   update->whichflag = 0;
   update->firststep = 0;
   update->laststep = nb_step;
   update->ntimestep = 0;
   update->beginstep = 0;
   update->endstep = update->laststep;
 
   if (flag_units) {
     force->ftm2v = code_unit_system.ft_m2v;
     // force->nktv2p = 1;
     force->mvv2e = code_unit_system.mvv2e;
   }
 
   if (this->domain->dimension != Dim)
     LM_FATAL("mismatch dimension: check lammps config file");
 
   this->acceleration_field.resize(atom->nlocal, 3);
   this->acceleration_field.setZero();
 
   if (std::string(this->atom->atom_style) != "atomic") {
     this->is_mass_per_particle = true;
     this->angular_acceleration_field.resize(atom->nlocal, 3);
     this->angular_acceleration_field.setZero();
   }
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> void DomainLammps<Dim>::initDOFs() {
   {
     UInt nall;
     if (force->newton)
       nall = atom->nlocal + atom->nghost;
     else
       nall = atom->nlocal;
 
     Real **x = atom->x;
     Real **x0 = atom->x0;
 
     for (UInt i = 0; i < nall; i++) {
       x0[i][0] = x[i][0];
       x0[i][1] = x[i][1];
       x0[i][2] = x[i][2];
     }
   }
   if (this->shouldRestart()) {
     this->restart_start = true;
     this->readXMLFile(this->restartfile);
 
     if (this->newGeomBox != invalidGeom &&
         current_step >= this->when_change_box) {
       this->changeBox();
       this->enforceCompatibility(); // perform migration after change box to
                                     // keep periodicity
       this->when_change_box = current_step + nb_step; // to avoid repetitions
     } else if (neighbor->decide()) {
       this->enforceCompatibility();
       this->resetBox();
     }
     this->updateGradient();
 
     LAMMPS_NS::LAMMPS::output->setup(false);
   }
   // conversion of masses from g/mol to Kg
   if (flag_units) {
     UnitsConverter uc;
     uc.setOutUnits(code_unit_system);
     if (strcmp(update->unit_style, "real") == 0)
       uc.setInUnits(real_unit_system);
     else if (strcmp(update->unit_style, "metal") == 0)
       uc.setInUnits(metal_unit_system);
     else if (strcmp(update->unit_style, "si") == 0)
       uc.setInUnits(si_unit_system);
     else
       LM_FATAL("not known lammps unit system:"
                << " things need to be done to integrate");
 
     uc.computeConversions();
 
     if (not is_mass_per_particle) {
       for (int i = 1; i <= atom->ntypes; i++)
         atom->mass[i] *= uc.getConversion<Mass>();
     } else {
       for (int i = 0; i <= atom->natoms; i++)
         atom->rmass[i] *= uc.getConversion<Mass>();
     }
   }
 
 #ifdef LIBMULTISCALE_TRACE_ATOM
   for (int i = 0; i < atom->nlocal; ++i) {
     Vector<3> X = {atom->x0[i][0], atom->x0[i][1], atom->x0[i][2]};
 
     SET_INTERNAL_TRACE_INDEX(X, i);
     IF_TRACED(X,
               "traced atom has been detected at position " << internal_index);
   }
   VIEW_ATOM(RefLammps<Dim>);
 #endif
 }
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim> Real DomainLammps<Dim>::potentialEnergy() {
   return update->minimize->pe_compute->compute_scalar();
 }
 
 /* -------------------------------------------------------------------------- */
 extern "C" {
 static UInt fargc = 1;
 static char prog[20] = "lammps";
 static char *___tmp = &prog[0];
 static char **fargv = &___tmp;
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> DomainLammps<Dim>::~DomainLammps() {}
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim>
 DomainLammps<Dim>::DomainLammps(const LMID &ID, CommGroup &group)
     : LMObject(ID), DomainAtomic<ContainerLammps<Dim>, Dim>(ID, group),
       LAMMPS(fargc, fargv, group.getMPIComm()), flag_units(true), triclinic(0) {
 
   newGeomBox = invalidGeom;
   when_change_box = 0;
   lattice = "";
   lattice_size = 0.;
   for (UInt i = 0; i < 6; ++i)
     replica[i] = 0;
   create_header_flag = false;
   restart_start = false;
 
   for (UInt i = 0; i < 3; ++i) {
     lattice_orient_x[i] = 0;
     lattice_orient_y[i] = 0;
     lattice_orient_z[i] = 0;
     lattice_origin[i] = 0.;
   }
   lattice_orient_x[0] = 1;
   lattice_orient_y[1] = 1;
   lattice_orient_z[2] = 1;
 
   lattice_spacing[0] = 1.0;
   lattice_spacing[1] = 1.0;
   lattice_spacing[2] = 1.0;
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> Cube DomainLammps<Dim>::getSubDomainBoundingBox() const {
   Cube c;
   c.setXmin(VectorView<3>(domain->sublo));
   c.setXmax(VectorView<3>(domain->subhi));
   return c;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim>
 void DomainLammps<Dim>::setSubDomainBoundingBox(const Cube &cube) {
 
   VectorView<3>(domain->sublo) = cube.getXmin();
   VectorView<3>(domain->subhi) = cube.getXmax();
 
   if (neighbor->decide())
     this->resetBox();
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> Cube DomainLammps<Dim>::getDomainBoundingBox() const {
 
   Cube c;
 
   if (!this->getCommGroup().amIinGroup())
     return c;
 
   c.setXmin(VectorView<3>(domain->boxlo));
   c.setXmax(VectorView<3>(domain->boxhi));
   return c;
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim>
 void DomainLammps<Dim>::setDomainBoundingBox(const Cube &cube) {
   VectorView<3>(domain->boxlo) = cube.getXmin();
   VectorView<3>(domain->boxhi) = cube.getXmax();
 
   if (neighbor->decide())
     this->resetBox();
 }
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim> void DomainLammps<Dim>::resetBox() {
   atom->setup();
   if (domain->triclinic)
     domain->x2lamda(atom->nlocal);
   domain->pbc();
   domain->reset_box();
   comm->setup();
   if (neighbor->style)
     neighbor->setup_bins();
   comm->exchange();
   if (atom->sortfreq > 0)
     atom->sort();
   comm->borders();
   if (domain->triclinic)
     domain->lamda2x(atom->nlocal + atom->nghost);
   neighbor->build(false);
 }
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim> void DomainLammps<Dim>::enforceCompatibility() {
   if (this->getCommGroup().amIinGroup()) {
     this->echangeAtomes();
     this->acceleration_field.resize(atom->nlocal, 3);
   }
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> void DomainLammps<Dim>::echangeAtomes() {
   VIEW_ATOM(RefLammps<Dim>);
   UInt nflag = neighbor->decide();
 
   // if (nflag == 0 && this->restart_start==false && this->update_pad == false)
   // {
   if (nflag == 0 && this->restart_start == false) {
     timer->stamp();
     comm->forward_comm();
     timer->stamp(LAMMPS_NS::Timer::COMM);
   } else {
     if (modify->n_pre_exchange)
       modify->pre_exchange();
     if (domain->triclinic)
       domain->x2lamda(atom->nlocal);
 
     domain->pbc();
     if (domain->box_change) {
       domain->reset_box();
       comm->setup();
       if (neighbor->style)
         neighbor->setup_bins();
     }
     timer->stamp();
 
     comm->exchange();
     comm->borders();
     if (domain->triclinic)
       domain->lamda2x(atom->nlocal + atom->nghost);
     timer->stamp(LAMMPS_NS::Timer::COMM);
     if (modify->n_pre_neighbor)
       modify->pre_neighbor();
     neighbor->build(0);
     timer->stamp(LAMMPS_NS::Timer::NEIGH);
     this->restart_start = false;
   }
 
   if (this->when_change_box == current_time) {
     // if one change box size, this part has to be fulfilled.
     if (modify->n_pre_exchange)
       modify->pre_exchange();
     if (domain->triclinic)
       domain->x2lamda(atom->nlocal);
     domain->pbc();
     if (domain->box_change) {
       domain->reset_box();
       comm->setup();
       if (neighbor->style)
         neighbor->setup_bins();
     }
     timer->stamp();
     comm->exchange();
     comm->borders();
     if (domain->triclinic)
       domain->lamda2x(atom->nlocal + atom->nghost);
     timer->stamp(LAMMPS_NS::Timer::COMM);
     if (modify->n_pre_neighbor)
       modify->pre_neighbor();
     neighbor->build(false);
     timer->stamp(LAMMPS_NS::Timer::NEIGH);
     this->restart_start = false;
   }
   VIEW_ATOM(RefLammps<Dim>);
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> bool DomainLammps<Dim>::isPeriodic(UInt dir) {
   switch (dir) {
   case 0:
     return this->domain->xperiodic;
     break;
   case 1:
     return this->domain->yperiodic;
     break;
   case 2:
     return this->domain->zperiodic;
     break;
   default:
     LM_FATAL("called on invalid-dimension");
     return false;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> void DomainLammps<Dim>::force_clear(UInt) {
   UInt i;
 
   // clear global force array
   // nall includes ghosts only if either newton flag is set
 
   UInt nall;
   if (force->newton)
     nall = atom->nlocal + atom->nghost;
   else
     nall = atom->nlocal;
 
   Real **f = atom->f;
   for (i = 0; i < nall; i++) {
     f[i][0] = 0.0;
     f[i][1] = 0.0;
     f[i][2] = 0.0;
   }
 
-  //   if (integ.torqueflag) {
-  //     Real **torque = atom->torque;
-  //     for (i = 0; i < nall; i++) {
-  //       torque[i][0] = 0.0;
-  //       torque[i][1] = 0.0;
-  //       torque[i][2] = 0.0;
-  //     }
-  //   }
+  if (lammps_main_object->atom->torque) {
+    Real **torque = atom->torque;
+    for (i = 0; i < nall; i++) {
+      torque[i][0] = 0.0;
+      torque[i][1] = 0.0;
+      torque[i][2] = 0.0;
+    }
+  }
 
   //   if (granflag) {
   //     Real **phia = atom->phia;
   //     for (i = 0; i < nall; i++) {
   //       phia[i][0] = 0.0;
   //       phia[i][1] = 0.0;
   //       phia[i][2] = 0.0;
   //     }
   //   }
 
   // clear pair force array if using it this timestep to compute virial
 
   //   if (vflag == 2 && pairflag) {
   //     if (atom->nmax > maxpair) {
   //       maxpair = atom->nmax;
   //       memory->destroy_2d_Real_array(f_pair);
   //       f_pair = memory->create_2d_Real_array(maxpair,3,"verlet:f_pair");
   //     }
   //     for (i = 0; i < nall; i++) {
   //       f_pair[i][0] = 0.0;
   //       f_pair[i][1] = 0.0;
   //       f_pair[i][2] = 0.0;
   //     }
   //   }
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> void DomainLammps<Dim>::setCurrentStep(UInt ts) {
   current_step = ts;
 
   // content of the Update::reset_timestep function
   update->eflag_global = update->vflag_global = ts;
   for (int i = 0; i < modify->ncompute; i++) {
     modify->compute[i]->invoked_scalar = -1;
     modify->compute[i]->invoked_vector = -1;
     modify->compute[i]->invoked_array = -1;
     modify->compute[i]->invoked_peratom = -1;
     modify->compute[i]->invoked_local = -1;
     // modify->compute[i]->eflag_global = -1;
   }
   update->ntimestep = ts;
 
   if (update->laststep < ts)
     update->laststep += ts;
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> void DomainLammps<Dim>::setTimeStep(Real ts) {
   update->dt = ts;
   for (int i = 0; i < modify->nfix; ++i)
     modify->fix[i]->init();
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> UInt DomainLammps<Dim>::getCurrentStep() {
   return update->ntimestep;
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> Real DomainLammps<Dim>::getTimeStep() { return update->dt; }
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim>
 int DomainLammps<Dim>::getOrientValue(UInt axis, UInt comp) {
   if (axis == 0) {
     return this->lattice_orient_x[comp];
   } else if (axis == 1) {
     return this->lattice_orient_y[comp];
   } else if (axis == 2) {
     return this->lattice_orient_z[comp];
   } else {
     LM_FATAL("wrong axis value: " << axis);
   }
 }
 /* -------------------------------------------------------------------------- */
 template <UInt Dim> Real DomainLammps<Dim>::getLatticeConstant() {
   return this->lattice_size;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim> void DomainLammps<Dim>::newmarkPredictor() {
   if (this->newGeomBox != invalidGeom &&
       current_step >= this->when_change_box) {
     this->changeBox();
   }
 
   VIEW_ATOM(RefLammps<Dim>);
   STARTTIMER("AtomStepPosition");
   this->update->ntimestep++;
   this->update->integrate->ev_set(this->update->ntimestep);
   this->modify->initial_integrate(0);
   this->echangeAtomes();
   STOPTIMER("AtomStepPosition");
   VIEW_ATOM(RefLammps<Dim>);
   this->changeRelease();
 }
 
 template <UInt Dim> void DomainLammps<Dim>::newmarkCorrector() {
   VIEW_ATOM(RefLammps<Dim>);
   DUMP("AtomStepVelocities x " << this->atom->x, DBG_INFO);
   STARTTIMER("AtomStepVelocities");
   this->modify->final_integrate();
   if (this->modify->n_end_of_step)
     this->modify->end_of_step();
 
   if (LAMMPS_NS::LAMMPS::output->next == this->update->ntimestep) {
     this->timer->stamp();
     LAMMPS_NS::LAMMPS::output->write(this->update->ntimestep);
     this->timer->stamp(LAMMPS_NS::Timer::OUTPUT);
   }
   STOPTIMER("AtomStepVelocities");
   VIEW_ATOM(RefLammps<Dim>);
   this->changeRelease();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim> void DomainLammps<Dim>::updateAcceleration() {
   if (!this->getCommGroup().amIinGroup())
     return;
 
   double *rmass = atom->rmass;
   double *mass = atom->mass;
   double *radius = atom->radius;
   double inertia = 0.4;
 
   int *type = atom->type;
   auto f = this->gradient();
 
   std::function<Real(UInt)> compute_mass_inverse;
 
   if (rmass)
     compute_mass_inverse = [&](UInt i) {
       return code_unit_system.ft_m2v / rmass[i];
     };
   else
     compute_mass_inverse = [&](UInt i) {
       return code_unit_system.ft_m2v / mass[type[i]];
     };
 
   for (UInt i = 0; i < UInt(lammps_main_object->atom->nlocal); ++i) {
     Real mass_inverse = compute_mass_inverse(i);
     this->acceleration_field(i, 0) = mass_inverse * f(i, 0);
     this->acceleration_field(i, 1) = mass_inverse * f(i, 1);
     this->acceleration_field(i, 2) = mass_inverse * f(i, 2);
   }
 
   if (std::string(this->atom->atom_style) != "atomic") {
     auto torque = this->getField(_torque);
     for (UInt i = 0; i < UInt(lammps_main_object->atom->nlocal); ++i) {
       Real mass_inverse = code_unit_system.ft_m2v / inertia /
                           (radius[i] * radius[i] * rmass[i]);
       this->angular_acceleration_field(i, 0) = mass_inverse * torque(i, 0);
       this->angular_acceleration_field(i, 1) = mass_inverse * torque(i, 1);
       this->angular_acceleration_field(i, 2) = mass_inverse * torque(i, 2);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt Dim> void DomainLammps<Dim>::updateGradient() {
   if (!this->getCommGroup().amIinGroup())
     return;
   VIEW_ATOM(RefLammps<Dim>);
   STARTTIMER("AtomStepForces");
   UInt eflag, vflag;
   eflag = this->update->integrate->eflag;
   vflag = this->update->integrate->vflag;
 
   this->force_clear(vflag);
   if (this->modify->n_pre_force) {
     this->modify->pre_force(vflag);
     timer->stamp(LAMMPS_NS::Timer::MODIFY);
   }
 
   if (this->force->pair && this->force->pair->compute_flag) {
     this->force->pair->compute(eflag, vflag);
     this->timer->stamp(LAMMPS_NS::Timer::PAIR);
   }
 
   if (this->atom->molecular) {
     if (this->force->bond)
       this->force->bond->compute(eflag, vflag);
     if (this->force->angle)
       this->force->angle->compute(eflag, vflag);
     if (this->force->dihedral)
       this->force->dihedral->compute(eflag, vflag);
     if (this->force->improper)
       this->force->improper->compute(eflag, vflag);
     this->timer->stamp(LAMMPS_NS::Timer::BOND);
   }
 
   if (this->force->kspace) {
     this->force->kspace->compute(eflag, vflag);
     this->timer->stamp(LAMMPS_NS::Timer::KSPACE);
   }
 
   if (modify->n_pre_reverse) {
     modify->pre_reverse(eflag, vflag);
     timer->stamp(LAMMPS_NS::Timer::MODIFY);
   }
 
   if (this->force->newton) {
     this->comm->reverse_comm();
     this->timer->stamp(LAMMPS_NS::Timer::COMM);
   }
 
   if (this->modify->n_post_force)
     this->modify->post_force(vflag);
 
   STOPTIMER("AtomStepForces");
   VIEW_ATOM(RefLammps<Dim>);
 
   this->changeRelease();
 }
 
 /* -------------------------------------------------------------------------- */
 /* LMDESC LAMMPS
    This plugin make the interface with lammps
 */
 /* LMHERITANCE domain_md */
 
 /* LM   EXAMPLE
    Section MultiScale RealUnits\\
    ...\\
    GEOMETRY myGeom CUBE BBOX -1 1 -1 1 -1 1\\
    MODEL LAMMPS md\\
    ...\\
    endSection\\ \\
 
    Section LAMMPS:md RealUnits \\
    LAMMPS_FILE lammps.in        \\
    DOMAIN_GEOMETRY 1            \\
    endSection
  */
 
 template <UInt Dim> void DomainLammps<Dim>::declareParams() {
 
   DomainAtomic<ContainerLammps<Dim>, Dim>::declareParams();
 
   /*LMKEYWORD LAMMPS_FILE
     This specifies the configuration file to be used for lammps.
     If the create header keyword is used then only the potential
     definition should be into that file. Otherwise, a complete
     lammps configuration file can be specified.
   */
 
   this->parseKeyword("LAMMPS_FILE", lammpsfile);
 
   /* LMKEYWORD CREATE_HEADER
      A header is automatically generated based on the information
      passed by keywords LATTICE, LATTICE\_SIZE, BOUNDARY and REPLICA
   */
   this->parseTag("CREATE_HEADER", create_header_flag, false);
 
   /* LMKEYWORD LATTICE
      Specifies the lattice to be used to generate the lammps header.
      Admissible values are lammps keywords
      (see \url{http://lammps.sandia.gov/doc/lattice.html})
   */
   this->parseKeyword("LATTICE", lattice, "");
 
   /* LMKEYWORD LATTICE_SIZE
      Lattice length parameter to generate the primitive lattice
   */
 
   this->parseKeyword("LATTICE_SIZE", lattice_size, 0.);
 
   /* LMKEYWORD LATTICE_ORIGIN
      Lattice length parameter to generate the primitive lattice
   */
 
   this->parseVectorKeyword("LATTICE_ORIGIN", 3, lattice_origin,
                            VEC_DEFAULTS(0., 0., 0.));
 
   /* LMKEYWORD LATTICE_ORIENTX
      Lattice length parameter to generate the primitive lattice
   */
 
   this->parseVectorKeyword("LATTICE_ORIENTX", 3, lattice_orient_x,
                            VEC_DEFAULTS(1, 0, 0));
 
   /* LMKEYWORD LATTICE_ORIENTY
      Lattice length parameter to generate the primitive lattice
   */
 
   this->parseVectorKeyword("LATTICE_ORIENTY", 3, lattice_orient_y,
                            VEC_DEFAULTS(0, 1, 0));
 
   /* LMKEYWORD LATTICE_ORIENTZ
      Lattice length parameter to generate the primitive lattice
   */
 
   this->parseVectorKeyword("LATTICE_ORIENTZ", 3, lattice_orient_z,
                            VEC_DEFAULTS(0, 0, 1));
 
   /* LMKEYWORD SPACING
      The minimal normalized lattice sizes
   */
 
   this->parseVectorKeyword("SPACING", 3, lattice_spacing,
                            VEC_DEFAULTS(1., 1., 1.));
 
   /* LMKEYWORD BOUNDARY
      Sequence of lammps code for boundary.
      (see \url{http://lammps.sandia.gov/doc/boundary.html})
   */
   this->parseVectorKeyword("BOUNDARY", 3, boundaries, VEC_DEFAULTS("", "", ""));
 
   /* LMKEYWORD REPLICA
      For the creation of atoms a region, specified as the number of
      replicas of the lattice must be provided, before calling the
 
      create atoms command.
   */
   this->parseVectorKeyword("REPLICA", 6, replica,
                            VEC_DEFAULTS(0, 0, 0, 0, 0, 0));
 
   /* LMKEYWORD CHANGE_DOMAIN_BOX
      This allows to specify live reshaping of the bounding box.
      Useful to enforce stresses states or dislocation displacement
      field with periodic boundary conditions.
   */
   this->parseKeyword("CHANGE_DOMAIN_BOX", newGeomBox, invalidGeom);
 
   /* LMKEYWORD CHANGE_DOMAIN_BOX_TIME
      Gives a timestep where the reshaping provided by CHANGE\_DOMAIN\_BOX
      should occur.
   */
   this->parseKeyword("CHANGE_DOMAIN_BOX_TIME", when_change_box, 0u);
 
   /* LMKEYWORD TRICLINIC
      This changes the orthogonal simulation box to triclinic box
   */
   this->parseKeyword("TRICLINIC", triclinic, 0);
 
   /* LMKEYWORD TILT
      This tils the simulation box by given amount value
      0: xy-direction; 1: xz-direction; 2: yz-direction
      (TILT 0 0.25: tilt simulation box in xy-direction by 0.25)
   */
   this->parseVectorKeyword("TILT", 2, tilt, VEC_DEFAULTS(0., 0.));
 
   /* LMKEYWORD UNITS_CONVERSION
      This is for debug purpose: force libmultiscale not to touch the
      units conversion factors of lammps.
   */
   this->parseTag("UNITS_CONVERSION", flag_units, true);
 }
 /* -------------------------------------------------------------------------- */
 
 template class DomainLammps<2>;
 template class DomainLammps<3>;
 __END_LIBMULTISCALE__