diff --git a/src/finish.cpp b/src/finish.cpp
index 1292cc7bf7..6a9b484444 100644
--- a/src/finish.cpp
+++ b/src/finish.cpp
@@ -43,7 +43,7 @@ void Finish::end(int flag)
 {
   int i,m,nneigh,nneighfull;
   int histo[10];
-  int loopflag,minflag,prdflag,timeflag,fftflag,histoflag,neighflag;
+  int loopflag,minflag,prdflag,tadflag,timeflag,fftflag,histoflag,neighflag;
   double time,tmp,ave,max,min;
   double time_loop,time_other;
   bigint natoms;
@@ -57,9 +57,10 @@ void Finish::end(int flag)
   // flag = 0 = just loop summary
   // flag = 1 = dynamics or minimization
   // flag = 2 = PRD
+  // flag = 3 = TAD
   
   loopflag = 1;
-  minflag = prdflag = timeflag = fftflag = histoflag = neighflag = 0;
+  minflag = prdflag = tadflag = timeflag = fftflag = histoflag = neighflag = 0;
 
   if (flag == 1) {
     if (update->whichflag == 2) minflag = 1;
@@ -69,6 +70,10 @@ void Finish::end(int flag)
   if (flag == 2) {
     prdflag = histoflag = neighflag = 1;
     
+  }
+  if (flag == 3) {
+    tadflag = histoflag = neighflag = 1;
+    
   }
 
   // loop stats
@@ -155,7 +160,6 @@ void Finish::end(int flag)
     }
   }
 
-
   // PRD stats using PAIR,BOND,KSPACE for dephase,dynamics,quench
 
   if (prdflag) {
@@ -216,6 +220,92 @@ void Finish::end(int flag)
     }
   }
 
+  // TAD stats using PAIR,BOND,KSPACE for neb,dynamics,quench
+
+  if (tadflag) {
+    if (me == 0) {
+      if (screen) fprintf(screen,"\n");
+      if (logfile) fprintf(logfile,"\n");
+    }
+
+    if (screen) fprintf(screen,"TAD stats:\n");
+    if (logfile) fprintf(logfile,"TAD stats:\n");
+
+    time = timer->array[TIME_PAIR];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"  NEB      time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"  NEB      time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+
+    time = timer->array[TIME_BOND];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"  Dynamics time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"  Dynamics time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+
+    time = timer->array[TIME_KSPACE];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"  Quench   time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"  Quench   time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+
+
+    time = timer->array[TIME_COMM];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"  Comm     time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"  Comm     time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+
+
+    time = timer->array[TIME_OUTPUT];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"  Output   time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"  Output   time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+
+    time = time_other;
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"  Other    time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"  Other    time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+  }
+
   // timing breakdowns
 
   if (timeflag) {
diff --git a/src/min_fire.cpp b/src/min_fire.cpp
index 93cffd6ca2..135e517427 100644
--- a/src/min_fire.cpp
+++ b/src/min_fire.cpp
@@ -106,6 +106,14 @@ int MinFire::iterate(int maxiter)
       vdotf += v[i][0]*f[i][0] + v[i][1]*f[i][1] + v[i][2]*f[i][2];
     MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,world);
 
+    // sum vdotf over replicas, if necessary
+    // this communicator would be invalid for multiprocess replicas
+
+    if (update->multireplica == 1) {
+      vdotf = vdotfall;
+      MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
+    }
+
     // if (v dot f) > 0:
     // v = (1-alpha) v + alpha |v| Fhat
     // |v| = length of v, Fhat = unit f
@@ -117,11 +125,28 @@ int MinFire::iterate(int maxiter)
       for (int i = 0; i < nlocal; i++)
 	vdotv += v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
       MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,world);
+
+      // sum vdotv over replicas, if necessary
+      // this communicator would be invalid for multiprocess replicas
+
+      if (update->multireplica == 1) {
+	vdotv = vdotvall;
+	MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
+      }
+
       fdotf = 0.0;
       for (int i = 0; i < nlocal; i++)
 	fdotf += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];
       MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,world);
 
+      // sum fdotf over replicas, if necessary
+      // this communicator would be invalid for multiprocess replicas
+      
+      if (update->multireplica == 1) {
+	fdotf = fdotfall;
+	MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
+      }
+
       scale1 = 1.0 - alpha;
       if (fdotfall == 0.0) scale2 = 0.0;
       else scale2 = alpha * sqrt(vdotvall/fdotfall);
@@ -162,6 +187,14 @@ int MinFire::iterate(int maxiter)
     }
     MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world);
 
+    // min dtv over replicas, if necessary
+    // this communicator would be invalid for multiprocess replicas
+
+    if (update->multireplica == 1) {
+      dtvone = dtv;
+      MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld);
+    }
+
     // Euler integration step
 
     double **x = atom->x;
diff --git a/src/min_quickmin.cpp b/src/min_quickmin.cpp
index dbdd65f675..89a3403d4e 100644
--- a/src/min_quickmin.cpp
+++ b/src/min_quickmin.cpp
@@ -102,6 +102,14 @@ int MinQuickMin::iterate(int maxiter)
       vdotf += v[i][0]*f[i][0] + v[i][1]*f[i][1] + v[i][2]*f[i][2];
     MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,world);
 
+    // sum vdotf over replicas, if necessary
+    // this communicator would be invalid for multiprocess replicas
+
+    if (update->multireplica == 1) {
+      vdotf = vdotfall;
+      MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
+    }
+
     if (vdotfall < 0.0) {
       last_negative = ntimestep;
       for (int i = 0; i < nlocal; i++)
@@ -112,6 +120,15 @@ int MinQuickMin::iterate(int maxiter)
       for (int i = 0; i < nlocal; i++)
 	fdotf += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];
       MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,world);
+
+      // sum fdotf over replicas, if necessary
+      // this communicator would be invalid for multiprocess replicas
+
+      if (update->multireplica == 1) {
+	fdotf = fdotfall;
+	MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
+      }
+
       if (fdotfall == 0.0) scale = 0.0;
       else scale = vdotfall/fdotfall;
       for (int i = 0; i < nlocal; i++) {
@@ -136,6 +153,14 @@ int MinQuickMin::iterate(int maxiter)
     }
     MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world);
 
+    // min dtv over replicas, if necessary
+    // this communicator would be invalid for multiprocess replicas
+
+    if (update->multireplica == 1) {
+      dtvone = dtv;
+      MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld);
+    }
+
     // Euler integration step
 
     double **x = atom->x;