From 69605295942797b0a9d1474f4a563b34ae845a5f Mon Sep 17 00:00:00 2001
From: robeme <robert.meissner@tuhh.de>
Date: Tue, 8 May 2018 10:02:10 +0200
Subject: [PATCH] set appropriate flags for computing a vector in constructor
 of fix restrain

---
 src/fix_restrain.cpp | 28 +++++++++++++++++++++++++---
 src/fix_restrain.h   |  2 ++
 2 files changed, 27 insertions(+), 3 deletions(-)

diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp
index 09cb8946e5..80775e0283 100644
--- a/src/fix_restrain.cpp
+++ b/src/fix_restrain.cpp
@@ -53,6 +53,9 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   global_freq = 1;
   extscalar = 1;
+  vector_flag = 1;
+  size_vector = 3;
+  extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
 
@@ -188,6 +191,10 @@ void FixRestrain::min_setup(int vflag)
 void FixRestrain::post_force(int vflag)
 {
   energy = 0.0;
+  
+  ebond = 0.0;
+  eangle = 0.0;
+  edihed = 0.0;
 
   for (int m = 0; m < nrestrain; m++)
     if (rstyle[m] == BOND) restrain_bond(m);
@@ -271,7 +278,8 @@ void FixRestrain::restrain_bond(int m)
   if (r > 0.0) fbond = -2.0*rk/r;
   else fbond = 0.0;
 
-  energy += rk*dr;
+  ebond += rk*dr;
+  energy += ebond;
 
   // apply force to each of 2 atoms
 
@@ -378,7 +386,8 @@ void FixRestrain::restrain_angle(int m)
   dtheta = acos(c) - target[m];
   tk = k * dtheta;
 
-  energy += tk*dtheta;
+  eangle += tk*dtheta;
+  energy += eangle;
 
   a = -2.0 * tk * s;
   a11 = a*c / rsq1;
@@ -558,7 +567,8 @@ void FixRestrain::restrain_dihedral(int m)
   df1 *= -mult[m];
   p += 1.0;
 
-  energy += k * p;
+  edihed += k * p;
+  energy += edihed;
 
   fg = vb1x*vb2xm + vb1y*vb2ym + vb1z*vb2zm;
   hg = vb3x*vb2xm + vb3y*vb2ym + vb3z*vb2zm;
@@ -635,3 +645,15 @@ double FixRestrain::compute_scalar()
   MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
   return energy_all;
 }
+
+/* ----------------------------------------------------------------------
+  return individual energy contributions
+------------------------------------------------------------------------- */
+
+double FixRestrain::compute_vector(int n)
+{
+  if (n == 0) return ebond;
+  if (n == 1) return eangle;
+  if (n == 2) return edihed;
+  return 0.0;
+}
diff --git a/src/fix_restrain.h b/src/fix_restrain.h
index 14699ed5af..c956e40982 100644
--- a/src/fix_restrain.h
+++ b/src/fix_restrain.h
@@ -36,6 +36,7 @@ class FixRestrain : public Fix {
   void post_force_respa(int, int, int);
   void min_post_force(int);
   double compute_scalar();
+  double compute_vector(int);
 
  private:
   int ilevel_respa;
@@ -46,6 +47,7 @@ class FixRestrain : public Fix {
   double *kstart,*kstop,*target;
   double *cos_target,*sin_target;
   double energy,energy_all;
+  double ebond,eangle,edihed;
 
   void restrain_bond(int);
   void restrain_angle(int);