added respa to fix_flow_gauss, not fully understood yet

src/USER-MISC/fix_flow_gauss.cpp

@@ -27,6 +27,7 @@
 #include "domain.h"
 #include "error.h"
 #include "citeme.h"
+#include "respa.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -63,6 +64,8 @@ FixFlowGauss::FixFlowGauss(LAMMPS *lmp, int narg, char **arg) :
   extvector = 1;
   size_vector = 3;
   global_freq = 1;    //data available every timestep
+  respa_level_support = 1; 
+  ilevel_respa = 0;  //default= innermost respa level
 
   dimension = domain->dimension;
 
@@ -109,9 +112,23 @@ int FixFlowGauss::setmask()
   int mask = 0;
   mask |= POST_FORCE;
   mask |= THERMO_ENERGY;
+  mask |= POST_FORCE_RESPA;
   return mask;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixFlowGauss::init()
+{
+  //if respa level specified by fix_modify, then override default here
+  //if specified level too high, set to max level
+  if (strstr(update->integrate_style,"respa")) {
+    int max_respa = ((Respa *) update->integrate)->nlevels-1;
+    if (respa_level >= 0) 
+      ilevel_respa = MIN(respa_level,max_respa);
+  }
+}
+
 /* ----------------------------------------------------------------------
    setup is called after the initial evaluation of forces before a run, so we
    must remove the total force here too
@@ -127,7 +144,13 @@ void FixFlowGauss::setup(int vflag)
   if (mTot <= 0.0)
     error->all(FLERR,"Invalid group mass in fix flow/gauss");
 
-  post_force(vflag);
+  if (strstr(update->integrate_style,"verlet"))
+    post_force(vflag);
+  else {
+    ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+    post_force_respa(vflag,ilevel_respa,0);
+    ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -201,6 +224,11 @@ void FixFlowGauss::post_force(int vflag)
 
 }
 
+void FixFlowGauss::post_force_respa(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == ilevel_respa) post_force(vflag);
+}
+
 /* ----------------------------------------------------------------------
    negative of work done by this fix
    This is only computed if requested, either with fix_modify energy yes, or with the energy keyword. Otherwise returns 0.

src/USER-MISC/fix_flow_gauss.h

@@ -30,10 +30,12 @@ FixStyle(flow/gauss,FixFlowGauss)
     public:
       FixFlowGauss(class LAMMPS *, int, char **);
       int setmask();
+      void init();
+      void setup(int);
+      void post_force(int);
+      void post_force_respa(int, int, int);
       double compute_scalar();
       double compute_vector(int n);
-      void post_force(int);
-      void setup(int);
 
     protected:
       int dimension;
@@ -44,6 +46,7 @@ FixStyle(flow/gauss,FixFlowGauss)
       double pe_tot;      //total added energy
       double dt;          //timestep
       bool workflag;      //if calculate work done by fix
+      int ilevel_respa;
 
     };
 

src/fix.cpp

@@ -139,7 +139,7 @@ void Fix::modify_params(int narg, char **arg)
     } else if (strcmp(arg[iarg],"respa") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
       if (!respa_level_support) error->all(FLERR,"Illegal fix_modify command");
-      int lvl = force->inumeric(FLERR,arg[1]);
+      int lvl = force->inumeric(FLERR,arg[iarg+1]);
       if (lvl < 0) error->all(FLERR,"Illegal fix_modify command");
       respa_level = lvl-1;
       iarg += 2;