From 8776b8103337cb30f5e7e509697fe85fb7077a69 Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Fri, 15 Sep 2017 00:46:18 -0400
Subject: [PATCH] add virial support to fix addforce

---
 src/fix_addforce.cpp | 44 ++++++++++++++++++++++++++++++++++++++------
 1 file changed, 38 insertions(+), 6 deletions(-)

diff --git a/src/fix_addforce.cpp b/src/fix_addforce.cpp
index 472d03f1f9..a05933e2fe 100644
--- a/src/fix_addforce.cpp
+++ b/src/fix_addforce.cpp
@@ -51,6 +51,7 @@ FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
   extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
+  virial_flag = 1;
 
   xstr = ystr = zstr = NULL;
 
@@ -237,10 +238,16 @@ void FixAddForce::post_force(int vflag)
   double **f = atom->f;
   int *mask = atom->mask;
   imageint *image = atom->image;
+  double v[6];
   int nlocal = atom->nlocal;
 
   if (update->ntimestep % nevery) return;
 
+  // energy and virial setup
+
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
+
   if (lmp->kokkos)
     atom->sync_modify(Host, (unsigned int) (F_MASK | MASK_MASK),
                       (unsigned int) F_MASK);
@@ -283,6 +290,14 @@ void FixAddForce::post_force(int vflag)
         f[i][0] += xvalue;
         f[i][1] += yvalue;
         f[i][2] += zvalue;
+        if (evflag) {
+          v[0] = xvalue * unwrap[0];
+          v[1] = yvalue * unwrap[1];
+          v[2] = zvalue * unwrap[2];
+          v[3] = xvalue * unwrap[1];
+          v[4] = xvalue * unwrap[2];
+          v[5] = yvalue * unwrap[2];
+          v_tally(i,v);
       }
 
   // variable force, wrap with clear/add
@@ -290,6 +305,7 @@ void FixAddForce::post_force(int vflag)
   // wrap with clear/add
 
   } else {
+    double unwrap[3];
 
     modify->clearstep_compute();
 
@@ -310,16 +326,32 @@ void FixAddForce::post_force(int vflag)
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
         if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;
+        domain->unmap(x[i],image[i],unwrap);
+        if (xstyle == ATOM) xvalue = sforce[i][0];
+        if (ystyle == ATOM) yvalue = sforce[i][1];
+        if (zstyle == ATOM) zvalue = sforce[i][2];
+
         if (estyle == ATOM) foriginal[0] += sforce[i][3];
+        else {
+          if (xstyle) foriginal[0] -= xvalue*unwrap[0];
+          if (ystyle) foriginal[0] -= yvalue*unwrap[1];
+          if (zstyle) foriginal[0] -= zvalue*unwrap[2];
+        }
         foriginal[1] += f[i][0];
         foriginal[2] += f[i][1];
         foriginal[3] += f[i][2];
-        if (xstyle == ATOM) f[i][0] += sforce[i][0];
-        else if (xstyle) f[i][0] += xvalue;
-        if (ystyle == ATOM) f[i][1] += sforce[i][1];
-        else if (ystyle) f[i][1] += yvalue;
-        if (zstyle == ATOM) f[i][2] += sforce[i][2];
-        else if (zstyle) f[i][2] += zvalue;
+        if (xstyle) f[i][0] += xvalue;
+        if (ystyle) f[i][1] += yvalue;
+        if (zstyle) f[i][2] += zvalue;
+        if (evflag) {
+          v[0] = xstyle ? fx*unwrap[0] : 0.0;
+          v[1] = ystyle ? fy*unwrap[1] : 0.0;
+          v[2] = zstyle ? fz*unwrap[2] : 0.0;
+          v[3] = xstyle ? fx*unwrap[1] : 0.0;
+          v[4] = xstyle ? fx*unwrap[2] : 0.0;
+          v[5] = ystyle ? fy*unwrap[2] : 0.0;
+          v_tally(i, v);
+        }
       }
   }
 }