From 3ab7be86648d5718a00dd833f2189d39b19152f0 Mon Sep 17 00:00:00 2001
From: kaymakme <kaymakme@msu.edu>
Date: Wed, 3 Mar 2021 20:36:06 -0500
Subject: [PATCH] QMMM interface related: ignore QM-MM bonds while forming bond
lists, update charge matrix so that the total charge constraint only applies
to QM atoms
---
sPuReMD/src/forces.c | 29 ++++++++++++++++++++++++++++-
1 file changed, 28 insertions(+), 1 deletion(-)
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 73179c2b..bb621573 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -498,6 +498,26 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
for ( i = 0; i < system->N_cm - 1; ++i )
{
+ #if defined(QMMM)
+ // total charge constraint on QM atoms
+ if ( system->atoms[i].qmmm_mask == TRUE )
+ {
+ H->j[*Htop] = i;
+ H->val[*Htop] = 1.0;
+
+ H_sp->j[*H_sp_top] = i;
+ H_sp->val[*H_sp_top] = 1.0;
+ }
+ else {
+ H->j[*Htop] = i;
+ H->val[*Htop] = 0.0;
+
+ H_sp->j[*H_sp_top] = i;
+ H_sp->val[*H_sp_top] = 0.0;
+ }
+ *Htop = *Htop + 1;
+ *H_sp_top = *H_sp_top + 1;
+ #else
H->j[*Htop] = i;
H->val[*Htop] = 1.0;
*Htop = *Htop + 1;
@@ -505,6 +525,7 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
H_sp->j[*H_sp_top] = i;
H_sp->val[*H_sp_top] = 1.0;
*H_sp_top = *H_sp_top + 1;
+ #endif
}
H->j[*Htop] = system->N_cm - 1;
@@ -856,7 +877,10 @@ static void Init_Forces( reax_system *system, control_params *control,
++H_sp_top;
}
}
-
+ #if defined(QMMM)
+ if ( system->atoms[i].qmmm_mask == TRUE && system->atoms[j].qmmm_mask == TRUE )
+ {
+ #endif
/* hydrogen bond lists */
if ( control->hbond_cut > 0.0
&& (ihb == H_ATOM || ihb == H_BONDING_ATOM)
@@ -1000,6 +1024,9 @@ static void Init_Forces( reax_system *system, control_params *control,
}
}
}
+#if defined(QMMM)
+ }
+#endif
}
/* diagonal entry */
--
GitLab