diff --git a/src/USER-MISC/README b/src/USER-MISC/README
index 4dfcb3bdac..2e433ef81d 100644
--- a/src/USER-MISC/README
+++ b/src/USER-MISC/README
@@ -43,6 +43,7 @@ fix pimd, Yuxing Peng (U Chicago), yuxing at uchicago.edu, 24 Nov 2014
 fix smd, Axel Kohlmeyer, akohlmey at gmail.com, 19 May 2008
 fix ti/rs, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013
 fix ti/spring, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013
+fix ttm/mod, Sergey Starikov and Vasily Pisarev (JIHT), pisarevvv at gmail.com, 2 Feb 2015
 improper_style cossq, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12
 improper_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12
 improper_style ring, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12
diff --git a/src/USER-MISC/fix_ttm_mod.cpp b/src/USER-MISC/fix_ttm_mod.cpp
new file mode 100644
index 0000000000..c72c2c97d0
--- /dev/null
+++ b/src/USER-MISC/fix_ttm_mod.cpp
@@ -0,0 +1,871 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing authors: (in addition to authors of original fix ttm)
+   Sergey Starikov (Joint Institute for High Temperatures of RAS)
+   Vasily Pisarev (Joint Institute for High Temperatures of RAS)
+------------------------------------------------------------------------- */
+
+#include "lmptype.h"
+#include "mpi.h"
+#include "math.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_ttm_mod.h"
+#include "atom.h"
+#include "force.h"
+#include "update.h"
+#include "domain.h"
+#include "region.h"
+#include "respa.h"
+#include "comm.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "error.h"
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+using namespace MathConst;
+
+#define MAXLINE 1024
+
+/* ---------------------------------------------------------------------- */
+
+FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 9) error->all(FLERR,"Illegal fix ttm/mod command");
+  vector_flag = 1;
+  size_vector = 2;
+  global_freq = 1;
+  extvector = 1;
+  nevery = 1;
+  restart_peratom = 1;
+  restart_global = 1;
+  seed = atoi(arg[3]);
+  nxnodes = atoi(arg[4]);
+  nynodes = atoi(arg[5]);
+  nznodes = atoi(arg[6]);
+  fpr = fopen(arg[7],"r");
+  if (fpr == NULL) {
+    char str[128];
+    sprintf(str,"Cannot open file %s",arg[7]);
+    error->one(FLERR,str);
+  }
+  fpr_2 = fopen(arg[8],"r");
+  if (fpr == NULL) {
+    char str[128];
+    sprintf(str,"Cannot open file %s",arg[8]);
+    error->one(FLERR,str);
+  }
+  nfileevery = atoi(arg[9]);
+  if (nfileevery) {
+    if (narg != 11) error->all(FLERR,"Illegal fix ttm/mod command");
+    MPI_Comm_rank(world,&me);
+    if (me == 0) {
+      fp = fopen(arg[10],"w");
+      if (fp == NULL) {
+        char str[128];
+        sprintf(str,"Cannot open fix ttm/mod file %s",arg[10]);
+        error->one(FLERR,str);
+      }
+    }
+  }
+  char linee[MAXLINE];
+  double tresh_d;
+  int tresh_i;
+  // C0 (metal)
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  esheat_0 = tresh_d;
+  // C1 (metal*10^3)
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  esheat_1 = tresh_d;
+  // C2 (metal*10^6)
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  esheat_2 = tresh_d;
+  // C3 (metal*10^9)
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  esheat_3 = tresh_d;
+  // C4 (metal*10^12)
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  esheat_4 = tresh_d;
+  // C_limit
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  C_limit = tresh_d;
+  //Temperature damping factor
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  T_damp = tresh_d;
+  // rho_e
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  electronic_density = tresh_d;
+  //thermal_diffusion
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  el_th_diff = tresh_d;
+  // gamma_p
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  gamma_p = tresh_d;
+  // gamma_s
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  gamma_s = tresh_d;
+  // v0
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  v_0 = tresh_d;
+  // average intensity of pulse (source of energy) (metal units)
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  intensity = tresh_d;
+  // coordinate of 1st surface in x-direction (in box units) - constant
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%d",&tresh_i);
+  surface_l = tresh_i;
+  // coordinate of 2nd surface in x-direction (in box units) - constant
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%d",&tresh_i);
+  surface_r = tresh_i;
+  // skin_layer = intensity is reduced (I=I0*exp[-x/skin_layer])
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%d",&tresh_i);
+  skin_layer = tresh_i;
+  // width of pulse (picoseconds)
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  width = tresh_d;
+  // factor of electronic pressure (PF) Pe = PF*Ce*Te
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  pres_factor = tresh_d;
+  // effective free path of electrons (angstrom)
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  free_path = tresh_d;
+  // ionic density (ions*angstrom^{-3})
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  ionic_density = tresh_d;
+  // if movsur = 0: surface is freezed
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%d",&tresh_i);
+  movsur = tresh_i;
+  // electron_temperature_min
+  fgets(linee,MAXLINE,fpr_2);
+  fgets(linee,MAXLINE,fpr_2);
+  sscanf(linee,"%lg",&tresh_d);
+  electron_temperature_min = tresh_d;
+  fclose(fpr_2);
+  //t_surface is determined by electronic temperature (not constant)
+  t_surface_l = surface_l;
+  mult_factor = intensity;
+  duration = 0.0;
+  surface_double = double(t_surface_l)*(domain->xprd/nxnodes);
+  if ((C_limit+esheat_0) < 0.0)
+    error->all(FLERR,"Fix ttm electronic_specific_heat must be >= 0.0");
+  if (electronic_density <= 0.0)
+    error->all(FLERR,"Fix ttm electronic_density must be > 0.0");
+  if (gamma_p < 0.0) error->all(FLERR,"Fix ttm gamma_p must be >= 0.0");
+  if (gamma_s < 0.0) error->all(FLERR,"Fix ttm gamma_s must be >= 0.0");
+  if (v_0 < 0.0) error->all(FLERR,"Fix ttm v_0 must be >= 0.0");
+  if (ionic_density <= 0.0) error->all(FLERR,"Fix ttm ionic_density must be > 0.0");
+  v_0_sq = v_0*v_0;
+  // error check
+  if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0)
+    error->all(FLERR,"Fix ttm number of nodes must be > 0");
+  if (seed <= 0) error->all(FLERR,"Invalid random number seed in fix ttm command");
+  if (surface_l < 0) error->all(FLERR,"Surface coordinates must be >= 0");
+  if (surface_l >= surface_r) error->all(FLERR, "Left surface coordinate must be less than right surface coordinate");
+  // initialize Marsaglia RNG with processor-unique seed
+  random = new RanMars(lmp,seed + comm->me);
+  // allocate per-type arrays for force prefactors
+  gfactor1 = new double[atom->ntypes+1];
+  gfactor2 = new double[atom->ntypes+1];
+  // allocate 3d grid variables
+  total_nnodes = nxnodes*nynodes*nznodes;
+  memory->create(nsum,nxnodes,nynodes,nznodes,"ttm/mod:nsum");
+  memory->create(nsum_all,nxnodes,nynodes,nznodes,"ttm/mod:nsum_all");
+  memory->create(T_initial_set,nxnodes,nynodes,nznodes,"ttm/mod:T_initial_set");
+  memory->create(sum_vsq,nxnodes,nynodes,nznodes,"ttm/mod:sum_vsq");
+  memory->create(sum_mass_vsq,nxnodes,nynodes,nznodes,"ttm/mod:sum_mass_vsq");
+  memory->create(sum_vsq_all,nxnodes,nynodes,nznodes,"ttm/mod:sum_vsq_all");
+  memory->create(sum_mass_vsq_all,nxnodes,nynodes,nznodes,
+                 "ttm/mod:sum_mass_vsq_all");
+  memory->create(T_electron_old,nxnodes,nynodes,nznodes,"ttm/mod:T_electron_old");
+  memory->create(T_electron,nxnodes,nynodes,nznodes,"ttm/mod:T_electron");
+  memory->create(net_energy_transfer,nxnodes,nynodes,nznodes,
+                 "ttm/mod:net_energy_transfer");
+  memory->create(net_energy_transfer_all,nxnodes,nynodes,nznodes,
+                 "ttm/mod:net_energy_transfer_all");
+  flangevin = NULL;
+  grow_arrays(atom->nmax);
+  // zero out the flangevin array
+  for (int i = 0; i < atom->nmax; i++) {
+    flangevin[i][0] = 0;
+    flangevin[i][1] = 0;
+    flangevin[i][2] = 0;
+  }
+  atom->add_callback(0);
+  atom->add_callback(1);
+  // set initial electron temperatures from user input file
+  if (me == 0) read_initial_electron_temperatures();
+  MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixTTMMod::~FixTTMMod()
+{
+  if (nfileevery && me == 0) fclose(fp);
+  delete random;
+  delete [] gfactor1;
+  delete [] gfactor2;
+  memory->destroy(nsum);
+  memory->destroy(nsum_all);
+  memory->destroy(T_initial_set);
+  memory->destroy(sum_vsq);
+  memory->destroy(sum_mass_vsq);
+  memory->destroy(sum_vsq_all);
+  memory->destroy(sum_mass_vsq_all);
+  memory->destroy(T_electron_first);
+  memory->destroy(T_electron_old);
+  memory->destroy(T_electron);
+  memory->destroy(flangevin);
+  memory->destroy(net_energy_transfer);
+  memory->destroy(net_energy_transfer_all);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixTTMMod::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  mask |= POST_FORCE_RESPA;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTMMod::init()
+{
+  if (domain->dimension == 2)
+    error->all(FLERR,"Cannot use fix ttm with 2d simulation");
+  if (domain->nonperiodic != 0)
+    error->all(FLERR,"Cannot use nonperiodic boundares with fix ttm");
+  if (domain->triclinic)
+    error->all(FLERR,"Cannot use fix ttm with triclinic box");
+  // set force prefactors
+  for (int i = 1; i <= atom->ntypes; i++) {
+    gfactor1[i] = - gamma_p / force->ftm2v;
+    gfactor2[i] =
+      sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
+  }
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        net_energy_transfer_all[ixnode][iynode][iznode] = 0;
+  if (strstr(update->integrate_style,"respa"))
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTMMod::setup(int vflag)
+{
+  if (strstr(update->integrate_style,"verlet"))
+    post_force_setup(vflag);
+  else {
+    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+    post_force_respa_setup(vflag,nlevels_respa-1,0);
+    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTMMod::post_force(int vflag)
+{
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  double dx = domain->xprd/nxnodes;
+  double dy = domain->yprd/nynodes;
+  double dz = domain->zprd/nynodes;
+  double gamma1,gamma2;
+  // apply damping and thermostat to all atoms in fix group
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
+      double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
+      double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
+      int ixnode = static_cast<int>(xscale*nxnodes);
+      int iynode = static_cast<int>(yscale*nynodes);
+      int iznode = static_cast<int>(zscale*nznodes);
+      while (ixnode > nxnodes-1) ixnode -= nxnodes;
+      while (iynode > nynodes-1) iynode -= nynodes;
+      while (iznode > nznodes-1) iznode -= nznodes;
+      while (ixnode < 0) ixnode += nxnodes;
+      while (iynode < 0) iynode += nynodes;
+      while (iznode < 0) iznode += nznodes;
+      if (T_electron[ixnode][iynode][iznode] < 0)
+        error->all(FLERR,"Electronic temperature dropped below zero");
+      double tsqrt = sqrt(T_electron[ixnode][iynode][iznode]);
+      gamma1 = gfactor1[type[i]];
+      double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+      if (vsq > v_0_sq) gamma1 *= (gamma_p + gamma_s)/gamma_p;
+      gamma2 = gfactor2[type[i]] * tsqrt;
+      if (ixnode >= surface_l){
+        if (ixnode < surface_r){
+          flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
+          flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
+          flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
+          double x_surf = dx*double(surface_l)+dx;
+          double x_at = x[i][0] - domain->boxlo[0];
+          int right_xnode = ixnode + 1;
+          int right_ynode = iynode + 1;
+          int right_znode = iznode + 1;
+          if (right_xnode == nxnodes) right_xnode = 0;
+          if (right_ynode == nynodes) right_ynode = 0;
+          if (right_znode == nznodes) right_znode = 0;
+          int left_xnode = ixnode - 1;
+          int left_ynode = iynode - 1;
+          int left_znode = iznode - 1;
+          if (left_xnode == -1) left_xnode = nxnodes - 1;
+          if (left_ynode == -1) left_ynode = nynodes - 1;
+          if (left_znode == -1) left_znode = nznodes - 1;
+          double T_i = T_electron[ixnode][iynode][iznode];
+          double T_il = T_electron[left_xnode][iynode][iznode];
+          double T_ir = T_electron[right_xnode][iynode][iznode];
+          double T_iu = T_electron[ixnode][right_ynode][iznode];
+          double T_if = T_electron[ixnode][iynode][right_znode];
+          double C_i = el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity;
+          double C_il = el_properties(T_electron[left_xnode][iynode][iznode]).el_heat_capacity;
+          double C_ir = el_properties(T_electron[right_xnode][iynode][iznode]).el_heat_capacity;
+          double C_iu = el_properties(T_electron[ixnode][right_ynode][iznode]).el_heat_capacity;
+          double C_if = el_properties(T_electron[ixnode][iynode][right_znode]).el_heat_capacity;
+          double diff_x = (x_at - x_surf)*(x_at - x_surf);
+          diff_x = pow(diff_x,0.5);
+          double len_factor = diff_x/(diff_x+free_path);
+          if (movsur == 1){
+            if (x_at >= x_surf){
+              flangevin[i][0] -= pres_factor/ionic_density*((C_ir*T_ir*free_path/(diff_x+free_path)/(diff_x+free_path)) +
+                                                            (len_factor/dx)*(C_ir*T_ir-C_i*T_i));
+              flangevin[i][1] -= pres_factor/ionic_density/dy*(C_iu*T_iu-C_i*T_i);
+              flangevin[i][2] -= pres_factor/ionic_density/dz*(C_if*T_if-C_i*T_i);
+            }
+          }
+          else{
+            flangevin[i][0] -= pres_factor/ionic_density/dx*(C_ir*T_ir-C_i*T_i);
+            flangevin[i][1] -= pres_factor/ionic_density/dy*(C_iu*T_iu-C_i*T_i);
+            flangevin[i][2] -= pres_factor/ionic_density/dz*(C_if*T_if-C_i*T_i);
+          }
+          f[i][0] += flangevin[i][0];
+          f[i][1] += flangevin[i][1];
+          f[i][2] += flangevin[i][2];
+        }
+      }
+      if (movsur == 1){
+        if (ixnode < surface_l){
+          t_surface_l = ixnode;
+        }
+      }
+    }
+  }
+  MPI_Allreduce(&t_surface_l,&surface_l,1,MPI_INT,MPI_MIN,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTMMod::post_force_setup(int vflag)
+{
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  // apply langevin forces that have been stored from previous run
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      f[i][0] += flangevin[i][0];
+      f[i][1] += flangevin[i][1];
+      f[i][2] += flangevin[i][2];
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTMMod::post_force_respa(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTMMod::post_force_respa_setup(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_force_setup(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTMMod::reset_dt()
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    gfactor2[i] =
+      sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
+}
+
+/* ----------------------------------------------------------------------
+   read in initial electron temperatures from a user-specified file
+   only called by proc 0
+------------------------------------------------------------------------- */
+
+void FixTTMMod::read_initial_electron_temperatures()
+{
+  char line[MAXLINE];
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        T_initial_set[ixnode][iynode][iznode] = 0;
+  // read initial electron temperature values from file
+  int ixnode,iynode,iznode;
+  double T_tmp;
+  while (1) {
+    if (fgets(line,MAXLINE,fpr) == NULL) break;
+    sscanf(line,"%d %d %d %lg",&ixnode,&iynode,&iznode,&T_tmp);
+    if (T_tmp < 0.0) error->one(FLERR,"Fix ttm electron temperatures must be >= 0.0");
+    T_electron[ixnode][iynode][iznode] = T_tmp;
+    T_initial_set[ixnode][iynode][iznode] = 1;
+  }
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        if (T_initial_set[ixnode][iynode][iznode] == 0)
+          error->one(FLERR,"Initial temperatures not all set in fix ttm");
+  // close file
+  fclose(fpr);
+}
+
+/* ---------------------------------------------------------------------- */
+
+el_heat_capacity_thermal_conductivity FixTTMMod::el_properties(double T_e)
+{
+  el_heat_capacity_thermal_conductivity properties;
+  double T_temp = T_e/1000.0, T_reduced = T_damp*T_temp;
+  double T2 = T_temp*T_temp;
+  double T3 = T2*T_temp;
+  double T4 = T3*T_temp;
+  double poly = esheat_0 + esheat_1*T_temp + esheat_2*T2 + esheat_3*T3 + esheat_4*T4;
+  properties.el_heat_capacity = electronic_density*(poly*exp(-T_reduced*T_reduced) + C_limit); // heat capacity
+  properties.el_thermal_conductivity = el_th_diff*properties.el_heat_capacity; // thermal conductivity
+  return properties;
+}
+double FixTTMMod::el_sp_heat_integral(double T_e)
+{
+  double T_temp = T_e/1000.0, T_reduced = T_damp*T_temp;
+  if (T_damp != 0)
+    return electronic_density*(MY_PIS*(3*esheat_4/pow(T_damp,5)+2*esheat_2/pow(T_damp,3)+4*esheat_0/T_damp)*erf(T_reduced)+
+                               4*esheat_3/pow(T_damp,4)+4*esheat_1/T_damp/T_damp-
+                               ((6*esheat_4*T_temp+4*esheat_3)/pow(T_damp,4)+
+                                (4*esheat_1+4*esheat_4*pow(T_temp,3)+4*esheat_3*T_temp*T_temp+4*esheat_2*T_temp)/T_damp/T_damp)*exp(-T_reduced*T_reduced))*125.0+electronic_density*C_limit*T_e;
+  else
+    return electronic_density*((esheat_0 + C_limit)*T_e + esheat_1*T_temp*T_e/2.0 + esheat_2*T_temp*T_temp*T_e/3.0 + esheat_3*pow(T_temp,3)*T_e/4.0 + esheat_4*pow(T_temp,4)*T_e/5.0);
+}
+void FixTTMMod::end_of_step()
+{
+  double **x = atom->x;
+  double **v = atom->v;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (movsur == 1){
+    for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+      for (int iynode = 0; iynode < nynodes; iynode++)
+        for (int iznode = 0; iznode < nznodes; iznode++){
+          double TTT = T_electron[ixnode][iynode][iznode];
+          if (TTT > 0){
+            if (ixnode < t_surface_l)
+              t_surface_l = ixnode;
+          }
+        }
+  }
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        net_energy_transfer[ixnode][iynode][iznode] = 0;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
+      double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
+      double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
+      int ixnode = static_cast<int>(xscale*nxnodes);
+      int iynode = static_cast<int>(yscale*nynodes);
+      int iznode = static_cast<int>(zscale*nznodes);
+      while (ixnode > nxnodes-1) ixnode -= nxnodes;
+      while (iynode > nynodes-1) iynode -= nynodes;
+      while (iznode > nznodes-1) iznode -= nznodes;
+      while (ixnode < 0) ixnode += nxnodes;
+      while (iynode < 0) iynode += nynodes;
+      while (iznode < 0) iznode += nznodes;
+      if (ixnode >= t_surface_l){
+        if (ixnode < surface_r)
+          net_energy_transfer[ixnode][iynode][iznode] +=
+            (flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
+             flangevin[i][2]*v[i][2]);
+      }
+    }
+  MPI_Allreduce(&net_energy_transfer[0][0][0],
+                &net_energy_transfer_all[0][0][0],
+                total_nnodes,MPI_DOUBLE,MPI_SUM,world);
+  double dx = domain->xprd/nxnodes;
+  double dy = domain->yprd/nynodes;
+  double dz = domain->zprd/nznodes;
+  double del_vol = dx*dy*dz;
+  double el_specific_heat = 0.0;
+  double el_thermal_conductivity = el_properties(electron_temperature_min).el_thermal_conductivity;
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        {
+          if (el_properties(T_electron[ixnode][iynode][iznode]).el_thermal_conductivity > el_thermal_conductivity)
+            el_thermal_conductivity = el_properties(T_electron[ixnode][iynode][iznode]).el_thermal_conductivity;
+          if (el_specific_heat > 0.0)
+            {
+              if ((T_electron[ixnode][iynode][iznode] > 0.0) && (el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity < el_specific_heat))
+                el_specific_heat = el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity;
+            }
+          else if (T_electron[ixnode][iynode][iznode] > 0.0) el_specific_heat = el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity;
+        }
+  // num_inner_timesteps = # of inner steps (thermal solves)
+  // required this MD step to maintain a stable explicit solve
+  int num_inner_timesteps = 1;
+  double inner_dt = update->dt;
+  double stability_criterion = 1.0 -
+    2.0*inner_dt/el_specific_heat *
+    (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
+  if (stability_criterion < 0.0) {
+    inner_dt = 0.25*el_specific_heat /
+      (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
+  }
+  num_inner_timesteps = static_cast<unsigned int>(update->dt/inner_dt) + 1;
+  inner_dt = update->dt/double(num_inner_timesteps);
+  if (num_inner_timesteps > 1000000)
+    error->warning(FLERR,"Too many inner timesteps in fix ttm",0);
+  for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;
+       ith_inner_timestep++) {
+    for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+      for (int iynode = 0; iynode < nynodes; iynode++)
+        for (int iznode = 0; iznode < nznodes; iznode++)
+          T_electron_old[ixnode][iynode][iznode] =
+            T_electron[ixnode][iynode][iznode];
+    // compute new electron T profile
+    duration = duration + inner_dt;
+    for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+      for (int iynode = 0; iynode < nynodes; iynode++)
+        for (int iznode = 0; iznode < nznodes; iznode++) {
+          int right_xnode = ixnode + 1;
+          int right_ynode = iynode + 1;
+          int right_znode = iznode + 1;
+          if (right_xnode == nxnodes) right_xnode = 0;
+          if (right_ynode == nynodes) right_ynode = 0;
+          if (right_znode == nznodes) right_znode = 0;
+          int left_xnode = ixnode - 1;
+          int left_ynode = iynode - 1;
+          int left_znode = iznode - 1;
+          if (left_xnode == -1) left_xnode = nxnodes - 1;
+          if (left_ynode == -1) left_ynode = nynodes - 1;
+          if (left_znode == -1) left_znode = nznodes - 1;
+          double skin_layer_d = double(skin_layer);
+          double ixnode_d = double(ixnode);
+          double surface_d = double(t_surface_l);
+          mult_factor = 0.0;
+          if (duration < width){
+            if (ixnode >= t_surface_l) mult_factor = (intensity/(dx*skin_layer_d))*exp((-1.0)*(ixnode_d - surface_d)/skin_layer_d);
+          }
+          if (ixnode < t_surface_l) net_energy_transfer_all[ixnode][iynode][iznode] = 0.0;
+          double cr_vac = 1;
+          if (T_electron_old[ixnode][iynode][iznode] == 0) cr_vac = 0;
+          double cr_v_l_x = 1;
+          if (T_electron_old[left_xnode][iynode][iznode] == 0) cr_v_l_x = 0;
+          double cr_v_r_x = 1;
+          if (T_electron_old[right_xnode][iynode][iznode] == 0) cr_v_r_x = 0;
+          double cr_v_l_y = 1;
+          if (T_electron_old[ixnode][left_ynode][iznode] == 0) cr_v_l_y = 0;
+          double cr_v_r_y = 1;
+          if (T_electron_old[ixnode][right_ynode][iznode] == 0) cr_v_r_y = 0;
+          double cr_v_l_z = 1;
+          if (T_electron_old[ixnode][iynode][left_znode] == 0) cr_v_l_z = 0;
+          double cr_v_r_z = 1;
+          if (T_electron_old[ixnode][iynode][right_znode] == 0) cr_v_r_z = 0;
+          if (cr_vac != 0) {
+            T_electron[ixnode][iynode][iznode] =
+              T_electron_old[ixnode][iynode][iznode] +
+              inner_dt/el_properties(T_electron_old[ixnode][iynode][iznode]).el_heat_capacity *
+              ((cr_v_r_x*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[right_xnode][iynode][iznode]/2.0).el_thermal_conductivity*
+                (T_electron_old[right_xnode][iynode][iznode]-T_electron_old[ixnode][iynode][iznode])/dx -
+                cr_v_l_x*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[left_xnode][iynode][iznode]/2.0).el_thermal_conductivity*
+                (T_electron_old[ixnode][iynode][iznode]-T_electron_old[left_xnode][iynode][iznode])/dx)/dx +
+               (cr_v_r_y*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][right_ynode][iznode]/2.0).el_thermal_conductivity*
+                (T_electron_old[ixnode][right_ynode][iznode]-T_electron_old[ixnode][iynode][iznode])/dy -
+                cr_v_l_y*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][left_ynode][iznode]/2.0).el_thermal_conductivity*
+                (T_electron_old[ixnode][iynode][iznode]-T_electron_old[ixnode][left_ynode][iznode])/dy)/dy +
+               (cr_v_r_z*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][iynode][right_znode]/2.0).el_thermal_conductivity*
+                (T_electron_old[ixnode][iynode][right_znode]-T_electron_old[ixnode][iynode][iznode])/dz -
+                cr_v_l_z*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][iynode][left_znode]/2.0).el_thermal_conductivity*
+                (T_electron_old[ixnode][iynode][iznode]-T_electron_old[ixnode][iynode][left_znode])/dz)/dz);
+            T_electron[ixnode][iynode][iznode]+=inner_dt/el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity*
+              (mult_factor -
+               net_energy_transfer_all[ixnode][iynode][iznode]/del_vol);
+          }
+          else T_electron[ixnode][iynode][iznode] =
+                 T_electron_old[ixnode][iynode][iznode];
+          if ((T_electron[ixnode][iynode][iznode] > 0.0) && (T_electron[ixnode][iynode][iznode] < electron_temperature_min))
+            T_electron[ixnode][iynode][iznode] = T_electron[ixnode][iynode][iznode] + 0.5*(electron_temperature_min - T_electron[ixnode][iynode][iznode]);
+        }
+  }
+  // output nodal temperatures for current timestep
+  if ((nfileevery) && !(update->ntimestep % nfileevery)) {
+    // compute atomic Ta for each grid point
+    for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+      for (int iynode = 0; iynode < nynodes; iynode++)
+        for (int iznode = 0; iznode < nznodes; iznode++) {
+          nsum[ixnode][iynode][iznode] = 0;
+          nsum_all[ixnode][iynode][iznode] = 0;
+          sum_vsq[ixnode][iynode][iznode] = 0.0;
+          sum_mass_vsq[ixnode][iynode][iznode] = 0.0;
+          sum_vsq_all[ixnode][iynode][iznode] = 0.0;
+          sum_mass_vsq_all[ixnode][iynode][iznode] = 0.0;
+        }
+    double massone;
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+        if (rmass) massone = rmass[i];
+        else massone = mass[type[i]];
+        double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
+        double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
+        double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
+        int ixnode = static_cast<int>(xscale*nxnodes);
+        int iynode = static_cast<int>(yscale*nynodes);
+        int iznode = static_cast<int>(zscale*nznodes);
+        while (ixnode > nxnodes-1) ixnode -= nxnodes;
+        while (iynode > nynodes-1) iynode -= nynodes;
+        while (iznode > nznodes-1) iznode -= nznodes;
+        while (ixnode < 0) ixnode += nxnodes;
+        while (iynode < 0) iynode += nynodes;
+        while (iznode < 0) iznode += nznodes;
+        double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+        nsum[ixnode][iynode][iznode] += 1;
+        sum_vsq[ixnode][iynode][iznode] += vsq;
+        sum_mass_vsq[ixnode][iynode][iznode] += massone*vsq;
+      }
+    MPI_Allreduce(&nsum[0][0][0],&nsum_all[0][0][0],total_nnodes,
+                  MPI_INT,MPI_SUM,world);
+    MPI_Allreduce(&sum_vsq[0][0][0],&sum_vsq_all[0][0][0],total_nnodes,
+                  MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(&sum_mass_vsq[0][0][0],&sum_mass_vsq_all[0][0][0],
+                  total_nnodes,MPI_DOUBLE,MPI_SUM,world);
+    MPI_Allreduce(&t_surface_l,&surface_l,
+                  1,MPI_INT,MPI_MIN,world);
+    if (me == 0) {
+      fprintf(fp,BIGINT_FORMAT,update->ntimestep);
+      double T_a;
+      for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+        for (int iynode = 0; iynode < nynodes; iynode++)
+          for (int iznode = 0; iznode < nznodes; iznode++) {
+            T_a = 0;
+            if (nsum_all[ixnode][iynode][iznode] > 0){
+              T_a = sum_mass_vsq_all[ixnode][iynode][iznode]/
+                (3.0*force->boltz*nsum_all[ixnode][iynode][iznode]/force->mvv2e);
+              if (movsur == 1){
+                if (T_electron[ixnode][iynode][iznode]==0.0) T_electron[ixnode][iynode][iznode] = electron_temperature_min;
+              }
+            }
+            fprintf(fp," %f",T_a);
+          }
+      fprintf(fp,"\t");
+      for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+        for (int iynode = 0; iynode < nynodes; iynode++)
+          for (int iznode = 0; iznode < nznodes; iznode++)
+            fprintf(fp,"%f ",T_electron[ixnode][iynode][iznode]);
+      fprintf(fp,"\n");
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of 3d grid
+------------------------------------------------------------------------- */
+
+double FixTTMMod::memory_usage()
+{
+  double bytes = 0.0;
+  bytes += 5*total_nnodes * sizeof(int);
+  bytes += 14*total_nnodes * sizeof(double);
+  return bytes;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixTTMMod::grow_arrays(int ngrow)
+{
+  memory->grow(flangevin,ngrow,3,"ttm/mod:flangevin");
+}
+
+/* ----------------------------------------------------------------------
+   return the energy of the electronic subsystem or the net_energy transfer
+   between the subsystems
+------------------------------------------------------------------------- */
+
+double FixTTMMod::compute_vector(int n)
+{
+  double e_energy = 0.0;
+  double transfer_energy = 0.0;
+  double dx = domain->xprd/nxnodes;
+  double dy = domain->yprd/nynodes;
+  double dz = domain->zprd/nznodes;
+  double del_vol = dx*dy*dz;
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++) {
+        e_energy += el_sp_heat_integral(T_electron[ixnode][iynode][iznode])*del_vol;
+        transfer_energy +=
+          net_energy_transfer_all[ixnode][iynode][iznode]*update->dt;
+      }
+  if (n == 0) return e_energy;
+  if (n == 1) return transfer_energy;
+  return 0.0;
+}
+
+/* ----------------------------------------------------------------------
+   pack entire state of Fix into one write
+------------------------------------------------------------------------- */
+
+void FixTTMMod::write_restart(FILE *fp)
+{
+  double *rlist;
+  memory->create(rlist,nxnodes*nynodes*nznodes+1,"ttm/mod:rlist");
+  int n = 0;
+  rlist[n++] = seed;
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        rlist[n++] = T_electron[ixnode][iynode][iznode];
+  if (comm->me == 0) {
+    int size = n * sizeof(double);
+    fwrite(&size,sizeof(int),1,fp);
+    fwrite(rlist,sizeof(double),n,fp);
+  }
+  memory->destroy(rlist);
+}
+
+/* ----------------------------------------------------------------------
+   use state info from restart file to restart the Fix
+------------------------------------------------------------------------- */
+
+void FixTTMMod::restart(char *buf)
+{
+  int n = 0;
+  double *rlist = (double *) buf;
+  // the seed must be changed from the initial seed
+  seed = static_cast<int> (0.5*rlist[n++]);
+  for (int ixnode = 0; ixnode < nxnodes; ixnode++)
+    for (int iynode = 0; iynode < nynodes; iynode++)
+      for (int iznode = 0; iznode < nznodes; iznode++)
+        T_electron[ixnode][iynode][iznode] = rlist[n++];
+  delete random;
+  random = new RanMars(lmp,seed+comm->me);
+}
+
+/* ----------------------------------------------------------------------
+   pack values in local atom-based arrays for restart file
+------------------------------------------------------------------------- */
+
+int FixTTMMod::pack_restart(int i, double *buf)
+{
+  buf[0] = 4;
+  buf[1] = flangevin[i][0];
+  buf[2] = flangevin[i][1];
+  buf[3] = flangevin[i][2];
+  return 4;
+}
+
+/* ----------------------------------------------------------------------
+   unpack values from atom->extra array to restart the fix
+------------------------------------------------------------------------- */
+
+void FixTTMMod::unpack_restart(int nlocal, int nth)
+{
+  double **extra = atom->extra;
+  // skip to Nth set of extra values
+  int m = 0;
+  for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
+  m++;
+  flangevin[nlocal][0] = extra[nlocal][m++];
+  flangevin[nlocal][1] = extra[nlocal][m++];
+  flangevin[nlocal][2] = extra[nlocal][m++];
+}
+
+/* ----------------------------------------------------------------------
+   maxsize of any atom's restart data
+------------------------------------------------------------------------- */
+
+int FixTTMMod::maxsize_restart()
+{
+  return 4;
+}
+
+/* ----------------------------------------------------------------------
+   size of atom nlocal's restart data
+------------------------------------------------------------------------- */
+
+int FixTTMMod::size_restart(int nlocal)
+{
+  return 4;
+}
diff --git a/src/USER-MISC/fix_ttm_mod.h b/src/USER-MISC/fix_ttm_mod.h
new file mode 100644
index 0000000000..11a68e4485
--- /dev/null
+++ b/src/USER-MISC/fix_ttm_mod.h
@@ -0,0 +1,88 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(ttm/mod,FixTTMMod)
+
+#else
+
+#ifndef LMP_FIX_TTM_MOD_H
+#define LMP_FIX_TTM_MOD_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+struct el_heat_capacity_thermal_conductivity {
+double el_heat_capacity;
+double el_thermal_conductivity;
+};
+
+class FixTTMMod : public Fix {
+ public:
+  FixTTMMod(class LAMMPS *, int, char **);
+  ~FixTTMMod();
+  int setmask();
+  void init();
+  void setup(int);
+  void post_force(int);
+  void post_force_respa(int, int, int);
+  void post_force_setup(int);
+  void post_force_respa_setup(int, int, int);
+  void end_of_step();
+  void reset_dt();
+  void write_restart(FILE *);
+  void restart(char *);
+  int pack_restart(int, double *);
+  void unpack_restart(int, int);
+  int size_restart(int);
+  int maxsize_restart();
+  double memory_usage();
+  void grow_arrays(int);
+  double compute_vector(int);
+
+ private:
+  int me;
+  int nfileevery;
+  int nlevels_respa;
+  int seed;
+  class RanMars *random;
+  FILE *fp,*fpr,*fpr_2;
+  int nxnodes,nynodes,nznodes,total_nnodes;
+  int ***nsum;
+  int ***nsum_all,***T_initial_set;
+  double *gfactor1,*gfactor2,*ratio;
+  double **flangevin;
+  double ***T_electron,***T_electron_old,***T_electron_first;
+  double ***sum_vsq,***sum_mass_vsq;
+  double ***sum_vsq_all,***sum_mass_vsq_all;
+  double ***net_energy_transfer,***net_energy_transfer_all;
+  double gamma_p,gamma_s,v_0,v_0_sq;
+  int skin_layer,surface_l,surface_r,t_surface_l,t_surface_r;
+  int movsur;
+  double esheat_0,esheat_1,esheat_2,esheat_3,esheat_4,C_limit,electronic_density;
+  double el_th_diff,T_damp;
+  double intensity,width,duration,surface_double;
+  double mult_factor,ttm_dt;
+  double pres_factor,free_path,ionic_density;
+  double electron_temperature_min;
+  el_heat_capacity_thermal_conductivity el_properties(double);
+  double el_sp_heat_integral(double);
+  void read_initial_electron_temperatures();
+};
+
+}
+
+#endif
+#endif