diff --git a/src/PHONON/dynamical_matrix.cpp b/src/PHONON/dynamical_matrix.cpp
index e13ab05b7f..d7c0ab17b0 100644
--- a/src/PHONON/dynamical_matrix.cpp
+++ b/src/PHONON/dynamical_matrix.cpp
@@ -49,7 +49,7 @@ enum{REGULAR,ESKM};
 
 DynamicalMatrix::DynamicalMatrix(LAMMPS *lmp) : Command(lmp), fp(nullptr)
 {
-    external_force_clear = 0;
+  external_force_clear = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -72,119 +72,118 @@ DynamicalMatrix::~DynamicalMatrix()
 
 void DynamicalMatrix::setup()
 {
-    // setup domain, communication and neighboring
-    // acquire ghosts
-    // build neighbor lists
-    if (triclinic) domain->x2lamda(atom->nlocal);
-    domain->pbc();
-    domain->reset_box();
-    if (neighbor->style) neighbor->setup_bins();
-    comm->exchange();
-    comm->borders();
-    if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
-    domain->image_check();
-    domain->box_too_small_check();
-    neighbor->build(1);
+  // setup domain, communication and neighboring
+  // acquire ghosts
+  // build neighbor lists
+  if (triclinic) domain->x2lamda(atom->nlocal);
+  domain->pbc();
+  domain->reset_box();
+  if (neighbor->style) neighbor->setup_bins();
+  comm->exchange();
+  comm->borders();
+  if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
+  domain->image_check();
+  domain->box_too_small_check();
+  neighbor->build(1);
 
-    // compute all forces
-    int ifix = modify->find_fix("package_omp");
-    if (ifix >= 0) external_force_clear = 1;
-    eflag=0;
-    vflag=0;
-    update_force();
+  // compute all forces
+  if (!modify->get_fix_by_id("package_omp")) external_force_clear = 1;
+  eflag=0;
+  vflag=0;
+  update_force();
 
-    if (pair_compute_flag) force->pair->compute(eflag,vflag);
-    else if (force->pair) force->pair->compute_dummy(eflag,vflag);
-    update->setupflag = 0;
+  if (pair_compute_flag) force->pair->compute(eflag,vflag);
+  else if (force->pair) force->pair->compute_dummy(eflag,vflag);
+  update->setupflag = 0;
 
-    //if all then skip communication groupmap population
-    if (gcount == atom->natoms)
-      for (bigint i=0; i<atom->natoms; i++)
-        groupmap[i] = i;
-    else
-      create_groupmap();
+  //if all then skip communication groupmap population
+  if (gcount == atom->natoms)
+    for (bigint i=0; i<atom->natoms; i++)
+      groupmap[i] = i;
+  else
+    create_groupmap();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DynamicalMatrix::command(int narg, char **arg)
 {
-    MPI_Comm_rank(world,&me);
+  MPI_Comm_rank(world,&me);
 
-    if (domain->box_exist == 0)
-        error->all(FLERR,"Dynamical_matrix command before simulation box is defined");
-    if (narg < 2) error->all(FLERR,"Illegal dynamical_matrix command");
+  if (domain->box_exist == 0)
+    error->all(FLERR,"Dynamical_matrix command before simulation box is defined");
+  if (narg < 2) error->all(FLERR,"Illegal dynamical_matrix command");
 
-    lmp->init();
+  lmp->init();
 
-    // orthogonal vs triclinic simulation box
+  // orthogonal vs triclinic simulation box
 
-    triclinic = domain->triclinic;
+  triclinic = domain->triclinic;
 
-    if (force->pair && force->pair->compute_flag) pair_compute_flag = 1;
-    else pair_compute_flag = 0;
-    if (force->kspace && force->kspace->compute_flag) kspace_compute_flag = 1;
-    else kspace_compute_flag = 0;
+  if (force->pair && force->pair->compute_flag) pair_compute_flag = 1;
+  else pair_compute_flag = 0;
+  if (force->kspace && force->kspace->compute_flag) kspace_compute_flag = 1;
+  else kspace_compute_flag = 0;
 
-    // group and style
+  // group and style
 
-    igroup = group->find(arg[0]);
-    if (igroup == -1) error->all(FLERR,"Could not find dynamical matrix group ID");
-    groupbit = group->bitmask[igroup];
-    gcount = group->count(igroup);
-    dynlen = (gcount)*3;
-    memory->create(groupmap,atom->natoms,"total_group_map:totalgm");
-    update->setupflag = 1;
+  igroup = group->find(arg[0]);
+  if (igroup == -1) error->all(FLERR,"Could not find dynamical matrix group ID");
+  groupbit = group->bitmask[igroup];
+  gcount = group->count(igroup);
+  dynlen = (gcount)*3;
+  memory->create(groupmap,atom->natoms,"total_group_map:totalgm");
+  update->setupflag = 1;
 
-    int style = -1;
-    if (strcmp(arg[1],"regular") == 0) style = REGULAR;
-    else if (strcmp(arg[1],"eskm") == 0) style = ESKM;
-    else error->all(FLERR,"Illegal Dynamical_matrix command");
-    del = utils::numeric(FLERR, arg[2],false,lmp);
+  int style = -1;
+  if (strcmp(arg[1],"regular") == 0) style = REGULAR;
+  else if (strcmp(arg[1],"eskm") == 0) style = ESKM;
+  else error->all(FLERR,"Illegal Dynamical_matrix command");
+  del = utils::numeric(FLERR, arg[2],false,lmp);
 
-    // set option defaults
+  // set option defaults
 
-    binaryflag = 0;
-    scaleflag = 0;
-    compressed = 0;
-    file_flag = 0;
-    file_opened = 0;
-    folded = 0;
-    conversion = 1;
+  binaryflag = 0;
+  scaleflag = 0;
+  compressed = 0;
+  file_flag = 0;
+  file_opened = 0;
+  folded = 0;
+  conversion = 1;
 
-    // read options from end of input line
-    if (style == REGULAR) options(narg-3,&arg[3]);
-    else if (style == ESKM) options(narg-3,&arg[3]);
-    else if (me == 0 && screen) fprintf(screen,"Illegal Dynamical Matrix command\n");
+  // read options from end of input line
+  if (style == REGULAR) options(narg-3,&arg[3]);
+  else if (style == ESKM) options(narg-3,&arg[3]);
+  else if (me == 0 && screen) fprintf(screen,"Illegal Dynamical Matrix command\n");
 
-    if (!folded) dynlenb = dynlen;
-    else dynlenb = (atom->natoms)*3;
+  if (!folded) dynlenb = dynlen;
+  else dynlenb = (atom->natoms)*3;
 
-    if (atom->map_style == Atom::MAP_NONE)
-        error->all(FLERR,"Dynamical_matrix command requires an atom map");
+  if (atom->map_style == Atom::MAP_NONE)
+    error->all(FLERR,"Dynamical_matrix command requires an atom map");
 
-    // move atoms by 3-vector or specified variable(s)
+  // move atoms by 3-vector or specified variable(s)
 
-    if (style == REGULAR) {
-        setup();
-        timer->init();
-        timer->barrier_start();
-        calculateMatrix();
-        timer->barrier_stop();
-    }
+  if (style == REGULAR) {
+    setup();
+    timer->init();
+    timer->barrier_start();
+    calculateMatrix();
+    timer->barrier_stop();
+  }
 
-    if (style == ESKM) {
-        setup();
-        convert_units(update->unit_style);
-        conversion = conv_energy/conv_distance/conv_mass;
-        timer->init();
-        timer->barrier_start();
-        calculateMatrix();
-        timer->barrier_stop();
-    }
+  if (style == ESKM) {
+    setup();
+    convert_units(update->unit_style);
+    conversion = conv_energy/conv_distance/conv_mass;
+    timer->init();
+    timer->barrier_start();
+    calculateMatrix();
+    timer->barrier_stop();
+  }
 
-    Finish finish(lmp);
-    finish.end(1);
+  Finish finish(lmp);
+  finish.end(1);
 }
 
 /* ----------------------------------------------------------------------
@@ -193,37 +192,37 @@ void DynamicalMatrix::command(int narg, char **arg)
 
 void DynamicalMatrix::options(int narg, char **arg)
 {
-    if (narg < 0) error->all(FLERR,"Illegal dynamical_matrix command");
-    int iarg = 0;
-    const char* filename = "dynmat.dyn";
+  if (narg < 0) error->all(FLERR,"Illegal dynamical_matrix command");
+  int iarg = 0;
+  const char* filename = "dynmat.dyn";
 
-    while (iarg < narg) {
-        if (strcmp(arg[iarg],"binary") == 0) {
-            if (iarg + 2 > narg) error->all(FLERR, "Illegal dynamical_matrix command");
-            if (strcmp(arg[iarg+1],"gzip") == 0) {
-                compressed = 1;
-            } else {
-                binaryflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
-            }
-            iarg += 2;
-      } else if (strcmp(arg[iarg],"file") == 0) {
-          if (iarg+2 > narg) error->all(FLERR, "Illegal dynamical_matrix command");
-          filename = arg[iarg + 1];
-          file_flag = 1;
-          iarg += 2;
-      } else if (strcmp(arg[iarg],"fold") == 0) {
-          if (iarg+2 > narg) error->all(FLERR, "Illegal dynamical_matrix command");
-          if (strcmp(arg[iarg+1],"yes") == 0) {
-              folded = 1;
-          } else if (strcmp(arg[iarg+1],"no") == 0) {
-              folded = 0;
-          } else error->all(FLERR,"Illegal input for dynamical_matrix fold option");
-          iarg += 2;
-      } else error->all(FLERR,"Illegal dynamical_matrix command");
-    }
-    if (file_flag == 1) {
-        openfile(filename);
-    }
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"binary") == 0) {
+      if (iarg + 2 > narg) error->all(FLERR, "Illegal dynamical_matrix command");
+      if (strcmp(arg[iarg+1],"gzip") == 0) {
+          compressed = 1;
+      } else {
+          binaryflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      }
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"file") == 0) {
+      if (iarg+2 > narg) error->all(FLERR, "Illegal dynamical_matrix command");
+      filename = arg[iarg + 1];
+      file_flag = 1;
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"fold") == 0) {
+      if (iarg+2 > narg) error->all(FLERR, "Illegal dynamical_matrix command");
+      if (strcmp(arg[iarg+1],"yes") == 0) {
+        folded = 1;
+      } else if (strcmp(arg[iarg+1],"no") == 0) {
+        folded = 0;
+      } else error->all(FLERR,"Illegal input for dynamical_matrix fold option");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal dynamical_matrix command");
+  }
+  if (file_flag == 1) {
+    openfile(filename);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -234,23 +233,23 @@ void DynamicalMatrix::options(int narg, char **arg)
 
 void DynamicalMatrix::openfile(const char *filename)
 {
-    // if file already opened, return
-    if (file_opened) return;
-    fp = nullptr;
+  // if file already opened, return
+  if (file_opened) return;
+  fp = nullptr;
 
-    if (me == 0) {
-        if (compressed) {
-            fp = platform::compressed_write(std::string(filename)+".gz");
-            if (!fp) error->one(FLERR,"Cannot open compressed file");
-        } else if (binaryflag) {
-            fp = fopen(filename,"wb");
-        } else {
-            fp = fopen(filename,"w");
-        }
-        if (!fp) error->one(FLERR,"Cannot open dynmat file: {}", utils::getsyserror());
+  if (me == 0) {
+    if (compressed) {
+      fp = platform::compressed_write(std::string(filename)+".gz");
+      if (!fp) error->one(FLERR,"Cannot open compressed file");
+    } else if (binaryflag) {
+      fp = fopen(filename,"wb");
+    } else {
+      fp = fopen(filename,"w");
     }
+    if (!fp) error->one(FLERR,"Cannot open dynmat file: {}", utils::getsyserror());
+  }
 
-    file_opened = 1;
+  file_opened = 1;
 }
 
 /* ----------------------------------------------------------------------
@@ -259,112 +258,112 @@ void DynamicalMatrix::openfile(const char *filename)
 
 void DynamicalMatrix::calculateMatrix()
 {
-    int local_idx; // local index
-    int local_jdx; // second local index
-    int nlocal = atom->nlocal;
-    bigint natoms = atom->natoms;
-    int *type = atom->type;
-    bigint *gm = groupmap;
-    double imass; // dynamical matrix element
-    double *m = atom->mass;
-    double **f = atom->f;
+  int local_idx; // local index
+  int local_jdx; // second local index
+  int nlocal = atom->nlocal;
+  bigint natoms = atom->natoms;
+  int *type = atom->type;
+  bigint *gm = groupmap;
+  double imass; // dynamical matrix element
+  double *m = atom->mass;
+  double **f = atom->f;
 
-    double **dynmat = new double*[3];
-    for (int i=0; i<3; i++)
-        dynmat[i] = new double[dynlenb];
+  double **dynmat = new double*[3];
+  for (int i=0; i<3; i++)
+    dynmat[i] = new double[dynlenb];
 
-    double **fdynmat = new double*[3];
-    for (int i=0; i<3; i++)
-        fdynmat[i] = new double[dynlenb];
+  double **fdynmat = new double*[3];
+  for (int i=0; i<3; i++)
+    fdynmat[i] = new double[dynlenb];
 
-    //initialize dynmat to all zeros
+  //initialize dynmat to all zeros
+  dynmat_clear(dynmat);
+
+  if (me == 0 && screen) {
+    fprintf(screen,"Calculating Dynamical Matrix ...\n");
+    fprintf(screen,"  Total # of atoms = " BIGINT_FORMAT "\n", natoms);
+    fprintf(screen,"  Atoms in group = " BIGINT_FORMAT "\n", gcount);
+    fprintf(screen,"  Total dynamical matrix elements = " BIGINT_FORMAT "\n", (dynlen*dynlen) );
+  }
+
+  // emit dynlen rows of dimalpha*dynlen*dimbeta elements
+
+  update->nsteps = 0;
+  int prog = 0;
+  for (bigint i=1; i<=natoms; i++) {
+    local_idx = atom->map(i);
+    if (gm[i-1] < 0)
+      continue;
+    for (int alpha=0; alpha<3; alpha++) {
+      displace_atom(local_idx, alpha, 1);
+      update_force();
+      for (bigint j=1; j<=natoms; j++) {
+        local_jdx = atom->map(j);
+        if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
+          && (gm[j-1] >= 0 || folded)){
+          if (folded) {
+            for (int beta=0; beta<3; beta++){
+              dynmat[alpha][(j-1)*3+beta] -= f[local_jdx][beta];
+            }
+          } else {
+            for (int beta=0; beta<3; beta++){
+              dynmat[alpha][gm[j-1]*3+beta] -= f[local_jdx][beta];
+            }
+          }
+        }
+      }
+      displace_atom(local_idx,alpha,-2);
+      update_force();
+      for (bigint j=1; j<=natoms; j++) {
+        local_jdx = atom->map(j);
+        if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
+          && (gm[j-1] >= 0 || folded)){
+          if (atom->rmass_flag == 1)
+            imass = sqrt(m[local_idx] * m[local_jdx]);
+          else
+            imass = sqrt(m[type[local_idx]] * m[type[local_jdx]]);
+          if (folded){
+            for (int beta=0; beta<3; beta++){
+              dynmat[alpha][(j-1)*3+beta] -= -f[local_jdx][beta];
+              dynmat[alpha][(j-1)*3+beta] /= (2 * del * imass);
+              dynmat[alpha][(j-1)*3+beta] *= conversion;
+            }
+          } else {
+            for (int beta=0; beta<3; beta++){
+              dynmat[alpha][gm[j-1]*3+beta] -= -f[local_jdx][beta];
+              dynmat[alpha][gm[j-1]*3+beta] /= (2 * del * imass);
+              dynmat[alpha][gm[j-1]*3+beta] *= conversion;
+            }
+          }
+        }
+      }
+      displace_atom(local_idx,alpha,1);
+    }
+    for (int k=0; k<3; k++)
+      MPI_Reduce(dynmat[k],fdynmat[k],dynlenb,MPI_DOUBLE,MPI_SUM,0,world);
+    if (me == 0)
+      writeMatrix(fdynmat);
     dynmat_clear(dynmat);
-
     if (me == 0 && screen) {
-        fprintf(screen,"Calculating Dynamical Matrix ...\n");
-        fprintf(screen,"  Total # of atoms = " BIGINT_FORMAT "\n", natoms);
-        fprintf(screen,"  Atoms in group = " BIGINT_FORMAT "\n", gcount);
-        fprintf(screen,"  Total dynamical matrix elements = " BIGINT_FORMAT "\n", (dynlen*dynlen) );
+      int p = 10 * gm[i-1] / gcount;
+      if (p > prog) {
+        prog = p;
+        fprintf(screen," %d%%",p*10);
+        fflush(screen);
+      }
     }
+  }
+  if (me == 0 && screen) fprintf(screen,"\n");
 
-    // emit dynlen rows of dimalpha*dynlen*dimbeta elements
+  for (int i=0; i < 3; i++)
+    delete [] dynmat[i];
+  delete [] dynmat;
 
-    update->nsteps = 0;
-    int prog = 0;
-    for (bigint i=1; i<=natoms; i++) {
-        local_idx = atom->map(i);
-        if (gm[i-1] < 0)
-            continue;
-        for (int alpha=0; alpha<3; alpha++) {
-            displace_atom(local_idx, alpha, 1);
-            update_force();
-            for (bigint j=1; j<=natoms; j++) {
-                local_jdx = atom->map(j);
-                if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
-                    && (gm[j-1] >= 0 || folded)){
-                    if (folded) {
-                        for (int beta=0; beta<3; beta++){
-                            dynmat[alpha][(j-1)*3+beta] -= f[local_jdx][beta];
-                        }
-                    } else {
-                        for (int beta=0; beta<3; beta++){
-                            dynmat[alpha][gm[j-1]*3+beta] -= f[local_jdx][beta];
-                        }
-                    }
-                }
-            }
-            displace_atom(local_idx,alpha,-2);
-            update_force();
-            for (bigint j=1; j<=natoms; j++) {
-                local_jdx = atom->map(j);
-                if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
-                    && (gm[j-1] >= 0 || folded)){
-                    if (atom->rmass_flag == 1)
-                        imass = sqrt(m[local_idx] * m[local_jdx]);
-                    else
-                        imass = sqrt(m[type[local_idx]] * m[type[local_jdx]]);
-                    if (folded){
-                        for (int beta=0; beta<3; beta++){
-                            dynmat[alpha][(j-1)*3+beta] -= -f[local_jdx][beta];
-                            dynmat[alpha][(j-1)*3+beta] /= (2 * del * imass);
-                            dynmat[alpha][(j-1)*3+beta] *= conversion;
-                        }
-                    } else {
-                        for (int beta=0; beta<3; beta++){
-                            dynmat[alpha][gm[j-1]*3+beta] -= -f[local_jdx][beta];
-                            dynmat[alpha][gm[j-1]*3+beta] /= (2 * del * imass);
-                            dynmat[alpha][gm[j-1]*3+beta] *= conversion;
-                        }
-                    }
-                }
-            }
-            displace_atom(local_idx,alpha,1);
-        }
-        for (int k=0; k<3; k++)
-            MPI_Reduce(dynmat[k],fdynmat[k],dynlenb,MPI_DOUBLE,MPI_SUM,0,world);
-        if (me == 0)
-            writeMatrix(fdynmat);
-        dynmat_clear(dynmat);
-        if (me == 0 && screen) {
-            int p = 10 * gm[i-1] / gcount;
-            if (p > prog) {
-                prog = p;
-                fprintf(screen," %d%%",p*10);
-                fflush(screen);
-            }
-        }
-    }
-    if (me == 0 && screen) fprintf(screen,"\n");
+  for (int i=0; i < 3; i++)
+    delete [] fdynmat[i];
+  delete [] fdynmat;
 
-    for (int i=0; i < 3; i++)
-        delete [] dynmat[i];
-    delete [] dynmat;
-
-    for (int i=0; i < 3; i++)
-        delete [] fdynmat[i];
-    delete [] fdynmat;
-
-    if (screen && me ==0 ) fprintf(screen,"Finished Calculating Dynamical Matrix\n");
+  if (screen && me ==0 ) fprintf(screen,"Finished Calculating Dynamical Matrix\n");
 }
 
 /* ----------------------------------------------------------------------
@@ -373,25 +372,25 @@ void DynamicalMatrix::calculateMatrix()
 
 void DynamicalMatrix::writeMatrix(double **dynmat)
 {
-    if (me != 0 || !fp)
-        return;
+  if (me != 0 || !fp)
+    return;
 
-    clearerr(fp);
-    if (binaryflag) {
-        for (int i=0; i<3; i++)
-            fwrite(dynmat[i], sizeof(double), dynlenb, fp);
-        if (ferror(fp))
-            error->one(FLERR, "Error writing to binary file");
-    } else {
-        for (int i = 0; i < 3; i++) {
-            for (bigint j = 0; j < dynlenb; j++) {
-                if ((j+1)%3==0) fprintf(fp, "%4.8f\n", dynmat[i][j]);
-                else fprintf(fp, "%4.8f ",dynmat[i][j]);
-            }
-        }
-        if (ferror(fp))
-            error->one(FLERR,"Error writing to file");
+  clearerr(fp);
+  if (binaryflag) {
+    for (int i=0; i<3; i++)
+      fwrite(dynmat[i], sizeof(double), dynlenb, fp);
+    if (ferror(fp))
+      error->one(FLERR, "Error writing to binary file");
+  } else {
+    for (int i = 0; i < 3; i++) {
+      for (bigint j = 0; j < dynlenb; j++) {
+        if ((j+1)%3==0) fprintf(fp, "%4.8f\n", dynmat[i][j]);
+        else fprintf(fp, "%4.8f ",dynmat[i][j]);
+      }
     }
+    if (ferror(fp))
+      error->one(FLERR,"Error writing to file");
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -400,18 +399,18 @@ void DynamicalMatrix::writeMatrix(double **dynmat)
 
 void DynamicalMatrix::displace_atom(int local_idx, int direction, int magnitude)
 {
-    if (local_idx < 0) return;
+  if (local_idx < 0) return;
 
-    double **x = atom->x;
-    int *sametag = atom->sametag;
-    int j = local_idx;
+  double **x = atom->x;
+  int *sametag = atom->sametag;
+  int j = local_idx;
 
-    x[local_idx][direction] += del*magnitude;
+  x[local_idx][direction] += del*magnitude;
 
-    while (sametag[j] >= 0) {
-        j = sametag[j];
-        x[j][direction] += del*magnitude;
-    }
+  while (sametag[j] >= 0) {
+    j = sametag[j];
+    x[j][direction] += del*magnitude;
+  }
 }
 
 
@@ -425,46 +424,46 @@ void DynamicalMatrix::displace_atom(int local_idx, int direction, int magnitude)
 
 void DynamicalMatrix::update_force()
 {
-    neighbor->decide(); // needed for intel potentials to work
-    force_clear();
-    int n_pre_force = modify->n_pre_force;
-    int n_pre_reverse = modify->n_pre_reverse;
-    int n_post_force = modify->n_post_force;
+  neighbor->decide(); // needed for intel potentials to work
+  force_clear();
+  int n_pre_force = modify->n_pre_force;
+  int n_pre_reverse = modify->n_pre_reverse;
+  int n_post_force = modify->n_post_force;
 
-    if (n_pre_force) {
-        modify->pre_force(vflag);
-        timer->stamp(Timer::MODIFY);
-    }
-
-    if (pair_compute_flag) {
-        force->pair->compute(eflag,vflag);
-        timer->stamp(Timer::PAIR);
-    }
-    if (atom->molecular != Atom::ATOMIC) {
-        if (force->bond) force->bond->compute(eflag,vflag);
-        if (force->angle) force->angle->compute(eflag,vflag);
-        if (force->dihedral) force->dihedral->compute(eflag,vflag);
-        if (force->improper) force->improper->compute(eflag,vflag);
-        timer->stamp(Timer::BOND);
-    }
-    if (kspace_compute_flag) {
-        force->kspace->compute(eflag,vflag);
-        timer->stamp(Timer::KSPACE);
-    }
-    if (n_pre_reverse) {
-        modify->pre_reverse(eflag,vflag);
-        timer->stamp(Timer::MODIFY);
-    }
-    if (force->newton) {
-        comm->reverse_comm();
-        timer->stamp(Timer::COMM);
-    }
-    // force modifications
-
-    if (n_post_force) modify->post_force(vflag);
+  if (n_pre_force) {
+    modify->pre_force(vflag);
     timer->stamp(Timer::MODIFY);
+  }
 
-    ++ update->nsteps;
+  if (pair_compute_flag) {
+    force->pair->compute(eflag,vflag);
+    timer->stamp(Timer::PAIR);
+  }
+  if (atom->molecular != Atom::ATOMIC) {
+    if (force->bond) force->bond->compute(eflag,vflag);
+    if (force->angle) force->angle->compute(eflag,vflag);
+    if (force->dihedral) force->dihedral->compute(eflag,vflag);
+    if (force->improper) force->improper->compute(eflag,vflag);
+    timer->stamp(Timer::BOND);
+  }
+  if (kspace_compute_flag) {
+    force->kspace->compute(eflag,vflag);
+    timer->stamp(Timer::KSPACE);
+  }
+  if (n_pre_reverse) {
+    modify->pre_reverse(eflag,vflag);
+    timer->stamp(Timer::MODIFY);
+  }
+  if (force->newton) {
+    comm->reverse_comm();
+    timer->stamp(Timer::COMM);
+  }
+  // force modifications
+
+  if (n_post_force) modify->post_force(vflag);
+  timer->stamp(Timer::MODIFY);
+
+  ++ update->nsteps;
 }
 
 /* ----------------------------------------------------------------------
@@ -474,17 +473,17 @@ void DynamicalMatrix::update_force()
 
 void DynamicalMatrix::force_clear()
 {
-    if (external_force_clear) return;
+  if (external_force_clear) return;
 
-    // clear global force array
-    // if either newton flag is set, also include ghosts
+  // clear global force array
+  // if either newton flag is set, also include ghosts
 
-    size_t nbytes = sizeof(double) * atom->nlocal;
-    if (force->newton) nbytes += sizeof(double) * atom->nghost;
+  size_t nbytes = sizeof(double) * atom->nlocal;
+  if (force->newton) nbytes += sizeof(double) * atom->nghost;
 
-    if (nbytes) {
-        memset(&atom->f[0][0],0,3*nbytes);
-    }
+  if (nbytes) {
+    memset(&atom->f[0][0],0,3*nbytes);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -493,67 +492,66 @@ void DynamicalMatrix::force_clear()
 
 void DynamicalMatrix::dynmat_clear(double **dynmat)
 {
+  size_t nbytes = sizeof(double) * dynlenb;
 
-    size_t nbytes = sizeof(double) * dynlenb;
-
-    if (nbytes) {
-        for (int i=0; i<3; i++)
-            memset(&dynmat[i][0],0,nbytes);
-    }
+  if (nbytes) {
+    for (int i=0; i<3; i++)
+      memset(&dynmat[i][0],0,nbytes);
+  }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void DynamicalMatrix::convert_units(const char *style)
 {
-    // physical constants from:
-    // https://physics.nist.gov/cuu/Constants/Table/allascii.txt
-    // using thermochemical calorie = 4.184 J
+  // physical constants from:
+  // https://physics.nist.gov/cuu/Constants/Table/allascii.txt
+  // using thermochemical calorie = 4.184 J
 
-    if (strcmp(style,"lj") == 0) {
-        error->all(FLERR,"Conversion Not Set");
-        //conversion = 1; // lj -> 10 J/mol
+  if (strcmp(style,"lj") == 0) {
+    error->all(FLERR,"Conversion Not Set");
+    //conversion = 1; // lj -> 10 J/mol
 
-    } else if (strcmp(style,"real") == 0) {
-        conv_energy = 418.4; // kcal/mol -> 10 J/mol
-        conv_mass = 1; // g/mol -> g/mol
-        conv_distance = 1; // angstrom -> angstrom
+  } else if (strcmp(style,"real") == 0) {
+    conv_energy = 418.4; // kcal/mol -> 10 J/mol
+    conv_mass = 1; // g/mol -> g/mol
+    conv_distance = 1; // angstrom -> angstrom
 
-    } else if (strcmp(style,"metal") == 0) {
-        conv_energy = 9648.5; // eV -> 10 J/mol
-        conv_mass = 1; // g/mol -> g/mol
-        conv_distance = 1; // angstrom -> angstrom
+  } else if (strcmp(style,"metal") == 0) {
+    conv_energy = 9648.5; // eV -> 10 J/mol
+    conv_mass = 1; // g/mol -> g/mol
+    conv_distance = 1; // angstrom -> angstrom
 
-    } else if (strcmp(style,"si") == 0) {
-        if (me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
-        conv_energy = 6.022E22; // J -> 10 J/mol
-        conv_mass = 6.022E26; // kg -> g/mol
-        conv_distance = 1E-10; // meter -> angstrom
+  } else if (strcmp(style,"si") == 0) {
+    if (me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
+    conv_energy = 6.022E22; // J -> 10 J/mol
+    conv_mass = 6.022E26; // kg -> g/mol
+    conv_distance = 1E-10; // meter -> angstrom
 
-    } else if (strcmp(style,"cgs") == 0) {
-        if (me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
-        conv_energy = 6.022E12; // Erg -> 10 J/mol
-        conv_mass = 6.022E23; // g -> g/mol
-        conv_distance = 1E-7; // centimeter -> angstrom
+  } else if (strcmp(style,"cgs") == 0) {
+    if (me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
+    conv_energy = 6.022E12; // Erg -> 10 J/mol
+    conv_mass = 6.022E23; // g -> g/mol
+    conv_distance = 1E-7; // centimeter -> angstrom
 
-    } else if (strcmp(style,"electron") == 0) {
-        conv_energy = 262550; // Hartree -> 10 J/mol
-        conv_mass = 1; // amu -> g/mol
-        conv_distance = 0.529177249; // bohr -> angstrom
+  } else if (strcmp(style,"electron") == 0) {
+    conv_energy = 262550; // Hartree -> 10 J/mol
+    conv_mass = 1; // amu -> g/mol
+    conv_distance = 0.529177249; // bohr -> angstrom
 
-    } else if (strcmp(style,"micro") == 0) {
-        if (me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
-        conv_energy = 6.022E10; // picogram-micrometer^2/microsecond^2 -> 10 J/mol
-        conv_mass = 6.022E11; // pg -> g/mol
-        conv_distance = 1E-4; // micrometer -> angstrom
+  } else if (strcmp(style,"micro") == 0) {
+    if (me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
+    conv_energy = 6.022E10; // picogram-micrometer^2/microsecond^2 -> 10 J/mol
+    conv_mass = 6.022E11; // pg -> g/mol
+    conv_distance = 1E-4; // micrometer -> angstrom
 
-    } else if (strcmp(style,"nano") == 0) {
-        if (me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
-        conv_energy = 6.022E4; // attogram-nanometer^2/nanosecond^2 -> 10 J/mol
-        conv_mass = 6.022E5; // ag -> g/mol
-        conv_distance = 0.1; // angstrom -> angstrom
+  } else if (strcmp(style,"nano") == 0) {
+    if (me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
+    conv_energy = 6.022E4; // attogram-nanometer^2/nanosecond^2 -> 10 J/mol
+    conv_mass = 6.022E5; // ag -> g/mol
+    conv_distance = 0.1; // angstrom -> angstrom
 
-    } else error->all(FLERR,"Units Type Conversion Not Found");
+  } else error->all(FLERR,"Units Type Conversion Not Found");
 
 }
 
@@ -561,66 +559,66 @@ void DynamicalMatrix::convert_units(const char *style)
 
 void DynamicalMatrix::create_groupmap()
 {
-    //Create a group map which maps atom order onto group
-    // groupmap[global atom index-1] = output column/row
+  //Create a group map which maps atom order onto group
+  // groupmap[global atom index-1] = output column/row
 
-    int local_idx; // local index
-    int gid = 0; //group index
-    int nlocal = atom->nlocal;
-    int *mask = atom->mask;
-    bigint natoms = atom->natoms;
-    int *recv = new int[comm->nprocs];
-    int *displs = new int[comm->nprocs];
-    bigint *temp_groupmap = new bigint[natoms];
+  int local_idx; // local index
+  int gid = 0; //group index
+  int nlocal = atom->nlocal;
+  int *mask = atom->mask;
+  bigint natoms = atom->natoms;
+  int *recv = new int[comm->nprocs];
+  int *displs = new int[comm->nprocs];
+  bigint *temp_groupmap = new bigint[natoms];
 
-    //find number of local atoms in the group (final_gid)
-    for (bigint i=1; i<=natoms; i++) {
-        local_idx = atom->map(i);
-        if ((local_idx >= 0) && (local_idx < nlocal) && mask[local_idx] & groupbit)
-            gid += 1; // gid at the end of loop is final_Gid
+  //find number of local atoms in the group (final_gid)
+  for (bigint i=1; i<=natoms; i++) {
+    local_idx = atom->map(i);
+    if ((local_idx >= 0) && (local_idx < nlocal) && mask[local_idx] & groupbit)
+      gid += 1; // gid at the end of loop is final_Gid
+  }
+  //create an array of length final_gid
+  bigint *sub_groupmap = new bigint[gid];
+
+  gid = 0;
+  //create a map between global atom id and group atom id for each proc
+  for (bigint i=1; i<=natoms; i++) {
+    local_idx = atom->map(i);
+    if ((local_idx >= 0) && (local_idx < nlocal) && mask[local_idx] & groupbit) {
+      sub_groupmap[gid] = i;
+      gid += 1;
     }
-    //create an array of length final_gid
-    bigint *sub_groupmap = new bigint[gid];
+  }
 
-    gid = 0;
-    //create a map between global atom id and group atom id for each proc
-    for (bigint i=1; i<=natoms; i++) {
-        local_idx = atom->map(i);
-        if ((local_idx >= 0) && (local_idx < nlocal) && mask[local_idx] & groupbit) {
-            sub_groupmap[gid] = i;
-            gid += 1;
-        }
-    }
+  //populate arrays for Allgatherv
+  for (int i=0; i < comm->nprocs; i++) {
+    recv[i] = 0;
+  }
+  recv[me] = gid;
+  MPI_Allreduce(recv,displs,comm->nprocs,MPI_INT,MPI_SUM,world);
+  for (int i=0; i<comm->nprocs; i++) {
+    recv[i]=displs[i];
+    if (i>0) displs[i] = displs[i-1]+recv[i-1];
+    else displs[i] = 0;
+  }
 
-    //populate arrays for Allgatherv
-    for (int i=0; i < comm->nprocs; i++) {
-        recv[i] = 0;
-    }
-    recv[me] = gid;
-    MPI_Allreduce(recv,displs,comm->nprocs,MPI_INT,MPI_SUM,world);
-    for (int i=0; i<comm->nprocs; i++) {
-        recv[i]=displs[i];
-        if (i>0) displs[i] = displs[i-1]+recv[i-1];
-        else displs[i] = 0;
-    }
+  //combine subgroup maps into total temporary groupmap
+  MPI_Allgatherv(sub_groupmap,gid,MPI_LMP_BIGINT,temp_groupmap,recv,displs,MPI_LMP_BIGINT,world);
+  std::sort(temp_groupmap,temp_groupmap+gcount);
 
-    //combine subgroup maps into total temporary groupmap
-    MPI_Allgatherv(sub_groupmap,gid,MPI_LMP_BIGINT,temp_groupmap,recv,displs,MPI_LMP_BIGINT,world);
-    std::sort(temp_groupmap,temp_groupmap+gcount);
+  //populate member groupmap based on temp groupmap
+  bigint j = 0;
+  for (bigint i=1; i<=natoms; i++) {
+    // flag groupmap contents that are in temp_groupmap
+    if (j < gcount && i == temp_groupmap[j])
+      groupmap[i-1] = j++;
+    else
+      groupmap[i-1] = -1;
+  }
 
-    //populate member groupmap based on temp groupmap
-    bigint j = 0;
-    for (bigint i=1; i<=natoms; i++) {
-        // flag groupmap contents that are in temp_groupmap
-        if (j < gcount && i == temp_groupmap[j])
-            groupmap[i-1] = j++;
-        else
-            groupmap[i-1] = -1;
-    }
-
-    //free that memory!
-    delete[] recv;
-    delete[] displs;
-    delete[] sub_groupmap;
-    delete[] temp_groupmap;
+  //free that memory!
+  delete[] recv;
+  delete[] displs;
+  delete[] sub_groupmap;
+  delete[] temp_groupmap;
 }