From 041df5c00e062a7eaedeccfd4564c1aaf0c99585 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Wed, 2 Jun 2021 14:49:47 -0400
Subject: [PATCH] sPuReMD: backport most recent QMMM changes for simulations
 with AmberMD.

---
 sPuReMD/src/charges.c   | 10 ++++++++--
 sPuReMD/src/forces.c    | 10 +++++++++-
 sPuReMD/src/neighbors.c | 16 ++++++++++++++++
 sPuReMD/src/spuremd.c   |  3 ---
 4 files changed, 33 insertions(+), 6 deletions(-)

diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index f3311c1a..998357a8 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -1518,6 +1518,9 @@ static void EE( reax_system * const system, control_params * const control,
         workspace->mask_qmmm[i] = system->atoms[i].qmmm_mask;
     }
     workspace->mask_qmmm[system->N_cm - 1] = 1;
+
+    /* Mask the b vector as well */
+    Vector_Mask_qmmm( workspace->b_s, workspace->mask_qmmm, system->N_cm );
 #endif
 
     switch ( control->cm_solver_type )
@@ -1646,8 +1649,11 @@ static void ACKS2( reax_system * const system, control_params * const control,
     {
         workspace->mask_qmmm[i] = system->atoms[i - system->N].qmmm_mask;
     }
-    workspace->mask_qmmm[2 * system->N_cm] = 1;
-    workspace->mask_qmmm[2 * system->N_cm + 1] = 1;
+    workspace->mask_qmmm[2 * system->N] = 1;
+    workspace->mask_qmmm[2 * system->N + 1] = 1;
+
+    /* Mask the b vector as well */
+    Vector_Mask_qmmm( workspace->b_s, workspace->mask_qmmm, system->N_cm );
 #endif
 
     switch ( control->cm_solver_type )
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 1b9cdc68..4b42eaa0 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -845,6 +845,11 @@ static void Init_Forces( reax_system *system, control_params *control,
             flag = FALSE;
             flag_sp = FALSE;
 
+#if defined(QMMM)
+            if ( system->atoms[i].qmmm_mask == TRUE
+                    || system->atoms[j].qmmm_mask == TRUE )
+            {
+#endif	
             /* check if reneighboring step --
              * atomic distances just computed via
              * Verlet list, so use current distances */
@@ -1076,11 +1081,14 @@ static void Init_Forces( reax_system *system, control_params *control,
 
                         Set_End_Index( j, btop_j + 1, bonds );
                     }
+#if defined(QMMM)
+                }
+#endif
                 }
-            }
 #if defined(QMMM)
             }
 #endif
+            }
         }
 
         /* diagonal entry */
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index a673c46d..8d6f5996 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -193,6 +193,11 @@ int Estimate_Num_Neighbors( reax_system * const system,
                             {
                                 atom2 = nbr_atoms[m];
 
+#if defined(QMMM)
+                                if ( system->atoms[atom1].qmmm_mask == TRUE
+                                        || system->atoms[atom2].qmmm_mask == TRUE )
+                                {
+#endif
                                 if ( atom1 >= atom2 )
                                 {
                                     count = Count_Far_Neighbors( system->atoms[atom1].x,
@@ -201,6 +206,9 @@ int Estimate_Num_Neighbors( reax_system * const system,
 
                                     num_far += count;
                                 }
+#if defined(QMMM)
+                                }
+#endif
                             }
                         }
 
@@ -298,6 +306,11 @@ void Generate_Neighbor_Lists( reax_system * const system,
                             {
                                 atom2 = nbr_atoms[m];
 
+#if defined(QMMM)
+                                if ( system->atoms[atom1].qmmm_mask == TRUE
+                                        || system->atoms[atom2].qmmm_mask == TRUE )
+                                {
+#endif
                                 if ( atom1 >= atom2 )
                                 {
                                     nbr_data = &far_nbrs->far_nbr_list[num_far];
@@ -308,6 +321,9 @@ void Generate_Neighbor_Lists( reax_system * const system,
 
                                     num_far += count;
                                 }
+#if defined(QMMM)
+                                }
+#endif
                             }
                         }
 
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 3664446c..7658e1a9 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -375,9 +375,6 @@ int simulate( const void * const handle )
         Reset( spmd_handle->system, spmd_handle->control, spmd_handle->data,
                 spmd_handle->workspace, spmd_handle->lists );
 
-        Generate_Neighbor_Lists( spmd_handle->system, spmd_handle->control, spmd_handle->data,
-                spmd_handle->workspace, spmd_handle->lists );
-
         Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
                 spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
                 spmd_handle->realloc );
-- 
GitLab