From 304156138347f16f77fcc15908ba9a9112cb822f Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Tue, 30 Mar 2021 07:36:53 -0400
Subject: [PATCH] add proper r-RESPA support for fix cmap

---
 doc/src/fix_cmap.rst      |  5 +++++
 src/MOLECULE/fix_cmap.cpp | 17 ++++++++++++-----
 src/MOLECULE/fix_cmap.h   |  2 +-
 3 files changed, 18 insertions(+), 6 deletions(-)

diff --git a/doc/src/fix_cmap.rst b/doc/src/fix_cmap.rst
index 892fd4ab41..e6e8d0c627 100644
--- a/doc/src/fix_cmap.rst
+++ b/doc/src/fix_cmap.rst
@@ -127,6 +127,11 @@ the :doc:`run <run>` command.
 The forces due to this fix are imposed during an energy minimization,
 invoked by the :doc:`minimize <minimize>` command.
 
+The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
+fix. This allows to set at which level of the :doc:`r-RESPA
+<run_style>` integrator the fix is adding its forces. Default is the
+outermost level.
+
 .. note::
 
    If you want the potential energy associated with the CMAP terms
diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp
index 6da259f6ca..1b58037eb7 100644
--- a/src/MOLECULE/fix_cmap.cpp
+++ b/src/MOLECULE/fix_cmap.cpp
@@ -81,6 +81,8 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) :
   extvector = 1;
   wd_header = 1;
   wd_section = 1;
+  respa_level_support = 1;
+  ilevel_respa = 0;
 
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
@@ -182,6 +184,11 @@ void FixCMAP::init()
   // define newton_bond here in case restart file was read (not data file)
 
   newton_bond = force->newton_bond;
+
+  if (utils::strmatch(update->integrate_style,"^respa")) {
+    ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+    if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+  }
 }
 
 /* --------------------------------------------------------------------- */
@@ -190,12 +197,12 @@ void FixCMAP::setup(int vflag)
 {
   pre_neighbor();
 
-  if (strstr(update->integrate_style,"verlet"))
+  if (utils::strmatch(update->integrate_style,"^verlet"))
     post_force(vflag);
   else {
-    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(vflag,nlevels_respa-1,0);
-    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+    ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
+    post_force_respa(vflag,ilevel_respa,0);
+    ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
   }
 }
 
@@ -596,7 +603,7 @@ void FixCMAP::post_force(int vflag)
 
 void FixCMAP::post_force_respa(int vflag, int ilevel, int /*iloop*/)
 {
-  if (ilevel == nlevels_respa-1) post_force(vflag);
+  if (ilevel == ilevel_respa) post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/MOLECULE/fix_cmap.h b/src/MOLECULE/fix_cmap.h
index 39d10bd77b..cd0c01a921 100644
--- a/src/MOLECULE/fix_cmap.h
+++ b/src/MOLECULE/fix_cmap.h
@@ -67,7 +67,7 @@ class FixCMAP : public Fix {
  private:
   int nprocs,me;
   int newton_bond,eflag_caller;
-  int ctype,nlevels_respa;
+  int ctype,ilevel_respa;
   int ncrosstermtypes,crossterm_per_atom,maxcrossterm;
   int ncrosstermlist;
   bigint ncmap;