diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 097c99c9e3ae8fd0bae863c3f34fa471ffa00f76..8ba64f9aff7c957a4d815fc7c538201e0f9ca98d 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -859,6 +859,10 @@ struct reax_system
int ffield_params_allocated;
/* number of local (non-periodic image) atoms for the current simulation */
int N;
+ /* number of local (non-periodic image) QM atoms for the current simulation in QMMM mode */
+ int N_qm;
+ /* number of local (non-periodic image) MM atoms for the current simulation in QMMM mode */
+ int N_mm;
/* max. number of local (non-periodic image) atoms across all simulations */
int N_max;
/* dimension of the N x N sparse charge method matrix H */
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 694dfc78c693f73b375a5288886a78bf642ef5ff..364b105b996ee01b0d973a1a1e63a2e361f9e213 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -227,6 +227,8 @@ void * setup_qmmm_( int qm_num_atoms, const int * const qm_types,
spmd_handle->callback = NULL;
spmd_handle->data->sim_id = 0;
+ spmd_handle->system->N_qm = qm_num_atoms;
+ spmd_handle->system->N_mm = mm_num_atoms;
spmd_handle->system->N = qm_num_atoms + mm_num_atoms;
PreAllocate_Space( spmd_handle->system, spmd_handle->control,
@@ -632,6 +634,8 @@ int reset_qmmm_( const void * const handle,
spmd_handle->realloc = FALSE;
spmd_handle->data->sim_id++;
+ spmd_handle->system->N_qm = qm_num_atoms;
+ spmd_handle->system->N_mm = mm_num_atoms;
spmd_handle->system->N = qm_num_atoms + mm_num_atoms;
PreAllocate_Space( spmd_handle->system, spmd_handle->control,
@@ -769,6 +773,181 @@ int reset( const void * const handle, const char * const geo_file,
}
+/* Getter for atom positions in QMMM mode
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * qm_pos_x: x-coordinate of QM atom positions, in Angstroms (allocated by caller)
+ * qm_pos_y: y-coordinate of QM atom positions, in Angstroms (allocated by caller)
+ * qm_pos_z: z-coordinate of QM atom positions, in Angstroms (allocated by caller)
+ * mm_pos_x: x-coordinate of MM atom positions, in Angstroms (allocated by caller)
+ * mm_pos_y: y-coordinate of MM atom positions, in Angstroms (allocated by caller)
+ * mm_pos_z: z-coordinate of MM atom positions, in Angstroms (allocated by caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int get_atom_positions_qmmm_( const void * const handle, double * const qm_pos_x,
+ double * const qm_pos_y, double * const qm_pos_z, double * const mm_pos_x,
+ double * const mm_pos_y, double * const mm_pos_z )
+{
+ int i, ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ for ( i = 0; i < spmd_handle->system->N_qm; ++i )
+ {
+ qm_pos_x[i] = spmd_handle->system->atoms[i].x[0];
+ qm_pos_y[i] = spmd_handle->system->atoms[i].x[1];
+ qm_pos_z[i] = spmd_handle->system->atoms[i].x[2];
+ }
+
+ for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
+ {
+ mm_pos_x[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].x[0];
+ mm_pos_y[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].x[1];
+ mm_pos_z[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].x[2];
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
+/* Getter for atom velocities in QMMM mode
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * qm_vel_x: x-coordinate of QM atom velocities, in Angstroms / ps (allocated by caller)
+ * qm_vel_y: y-coordinate of QM atom velocities, in Angstroms / ps (allocated by caller)
+ * qm_vel_z: z-coordinate of QM atom velocities, in Angstroms / ps (allocated by caller)
+ * mm_vel_x: x-coordinate of MM atom velocities, in Angstroms / ps (allocated by caller)
+ * mm_vel_y: y-coordinate of MM atom velocities, in Angstroms / ps (allocated by caller)
+ * mm_vel_z: z-coordinate of MM atom velocities, in Angstroms / ps (allocated by caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int get_atom_velocities_qmmm_( const void * const handle, double * const qm_vel_x,
+ double * const qm_vel_y, double * const qm_vel_z, double * const mm_vel_x,
+ double * const mm_vel_y, double * const mm_vel_z )
+{
+ int i, ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ for ( i = 0; i < spmd_handle->system->N_qm; ++i )
+ {
+ qm_vel_x[i] = spmd_handle->system->atoms[i].v[0];
+ qm_vel_y[i] = spmd_handle->system->atoms[i].v[1];
+ qm_vel_z[i] = spmd_handle->system->atoms[i].v[2];
+ }
+
+ for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
+ {
+ mm_vel_x[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].v[0];
+ mm_vel_y[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].v[1];
+ mm_vel_z[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].v[2];
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
+/* Getter for atom forces in QMMM mode
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * qm_f_x: x-coordinate of QM atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
+ * qm_f_y: y-coordinate of QM atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
+ * qm_f_z: z-coordinate of QM atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
+ * mm_f_x: x-coordinate of MM atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
+ * mm_f_y: y-coordinate of MM atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
+ * mm_f_z: z-coordinate of MM atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int get_atom_forces_qmmm_( const void * const handle, double * const qm_f_x,
+ double * const qm_f_y, double * const qm_f_z, double * const mm_f_x,
+ double * const mm_f_y, double * const mm_f_z )
+{
+ int i, ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ for ( i = 0; i < spmd_handle->system->N_qm; ++i )
+ {
+ qm_f_x[i] = spmd_handle->system->atoms[i].f[0];
+ qm_f_y[i] = spmd_handle->system->atoms[i].f[1];
+ qm_f_z[i] = spmd_handle->system->atoms[i].f[2];
+ }
+
+ for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
+ {
+ mm_f_x[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].f[0];
+ mm_f_y[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].f[1];
+ mm_f_z[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].f[2];
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
+/* Getter for atom charges in QMMM mode
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * qm_q: QM atom charges, in Coulombs (allocated by caller)
+ * mm_q: MM atom charges, in Coulombs (allocated by caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int get_atom_charges_qmmm_( const void * const handle, double * const qm_q,
+ double * const mm_q )
+{
+ int i, ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ for ( i = 0; i < spmd_handle->system->N_qm; ++i )
+ {
+ qm_q[i] = spmd_handle->system->atoms[i].q;
+ }
+
+ for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
+ {
+ mm_q[i - spmd_handle->system->N_qm] = spmd_handle->system->atoms[i].q;
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
/* Getter for atom positions
*
* handle: pointer to wrapper struct with top-level data structures
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index 3ee91704cbc2d00db5325f2e33d75b806ab46a57..74803b5afa7192f6fdcaddec1b182a4fc6386769 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -62,6 +62,20 @@ int reset_qmmm_( const void * const, int,
int reset( const void * const, const char * const,
const char * const, const char * const );
+int get_atom_positions_qmmm_( const void * const, double * const,
+ double * const, double * const, double * const,
+ double * const, double * const );
+
+int get_atom_velocities_qmmm_( const void * const, double * const,
+ double * const, double * const, double * const,
+ double * const, double * const );
+
+int get_atom_forces_qmmm_( const void * const, double * const,
+ double * const, double * const, double * const,
+ double * const, double * const );
+
+int get_atom_charges_qmmm_( const void * const, double * const, double * const );
+
int get_atom_positions( const void * const, double * const,
double * const, double * const );
diff --git a/tools/driver_qmmm.py b/tools/driver_qmmm.py
index ec315724f1c881ee357e3e7aa9ea74e5b545fe42..dfa80de68b7ec499cd119dc18894bc799da37e90 100644
--- a/tools/driver_qmmm.py
+++ b/tools/driver_qmmm.py
@@ -393,19 +393,21 @@ if __name__ == '__main__':
set_control_parameter.argtypes = [c_void_p, c_char_p, POINTER(c_char_p)]
set_control_parameter.restype = c_int
- get_atom_positions = lib.get_atom_positions
- get_atom_positions.argtypes = [c_void_p, POINTER(c_double),
- POINTER(c_double), POINTER(c_double)]
- get_atom_positions.restype = c_int
+ get_atom_positions_qmmm = lib.get_atom_positions_qmmm_
+ get_atom_positions_qmmm.argtypes = [c_void_p,
+ POINTER(c_double), POINTER(c_double), POINTER(c_double),
+ POINTER(c_double), POINTER(c_double), POINTER(c_double)]
+ get_atom_positions_qmmm.restype = c_int
- get_atom_forces = lib.get_atom_forces
- get_atom_forces.argtypes = [c_void_p, POINTER(c_double),
- POINTER(c_double), POINTER(c_double)]
- get_atom_forces.restype = c_int
+ get_atom_forces_qmmm = lib.get_atom_forces_qmmm_
+ get_atom_forces_qmmm.argtypes = [c_void_p,
+ POINTER(c_double), POINTER(c_double), POINTER(c_double),
+ POINTER(c_double), POINTER(c_double), POINTER(c_double)]
+ get_atom_forces_qmmm.restype = c_int
- get_atom_charges = lib.get_atom_charges
- get_atom_charges.argtypes = [c_void_p, POINTER(c_double)]
- get_atom_charges.restype = c_int
+ get_atom_charges_qmmm = lib.get_atom_charges_qmmm_
+ get_atom_charges_qmmm.argtypes = [c_void_p, POINTER(c_double), POINTER(c_double)]
+ get_atom_charges_qmmm.restype = c_int
def get_simulation_step_results(num_atoms, atoms, data):
print("{0:24.15f} {1:24.15f} {2:24.15f}".format(
@@ -450,27 +452,35 @@ if __name__ == '__main__':
if ret != 0:
print("[ERROR] simulate returned {0}".format(ret))
- p_x = (c_double * num_atoms)()
- p_y = (c_double * num_atoms)()
- p_z = (c_double * num_atoms)()
- ret = get_atom_positions(handle, p_x, p_y, p_z)
+ qm_p_x = (c_double * num_qm_atoms)()
+ qm_p_y = (c_double * num_qm_atoms)()
+ qm_p_z = (c_double * num_qm_atoms)()
+ mm_p_x = (c_double * num_mm_atoms)()
+ mm_p_y = (c_double * num_mm_atoms)()
+ mm_p_z = (c_double * num_mm_atoms)()
+ ret = get_atom_positions_qmmm(handle, qm_p_x, qm_p_y, qm_p_z,
+ mm_p_x, mm_p_y, mm_p_z)
if ret != 0:
- print("[ERROR] get_atom_positions returned {0}".format(ret))
+ print("[ERROR] get_atom_positions_qmmm returned {0}".format(ret))
- f_x = (c_double * num_atoms)()
- f_y = (c_double * num_atoms)()
- f_z = (c_double * num_atoms)()
- ret = get_atom_forces(handle, f_x, f_y, f_z)
+ qm_f_x = (c_double * num_qm_atoms)()
+ qm_f_y = (c_double * num_qm_atoms)()
+ qm_f_z = (c_double * num_qm_atoms)()
+ mm_f_x = (c_double * num_mm_atoms)()
+ mm_f_y = (c_double * num_mm_atoms)()
+ mm_f_z = (c_double * num_mm_atoms)()
+ ret = get_atom_forces_qmmm(handle, qm_f_x, qm_f_y, qm_f_z, mm_f_x, mm_f_y, mm_f_z)
if ret != 0:
- print("[ERROR] get_atom_forces returned {0}".format(ret))
+ print("[ERROR] get_atom_forces_qmmm returned {0}".format(ret))
- q = (c_double * num_atoms)()
- ret = get_atom_charges(handle, q)
+ qm_q = (c_double * num_qm_atoms)()
+ mm_q = (c_double * num_mm_atoms)()
+ ret = get_atom_charges_qmmm(handle, qm_q, mm_q)
if ret != 0:
- print("[ERROR] get_atom_charges returned {0}".format(ret))
+ print("[ERROR] get_atom_charges_qmmm returned {0}".format(ret))
# silica
sim_box = (c_double * 6)(36.477, 50.174, 52.110, 90.0, 90.0, 90.0)
@@ -499,26 +509,34 @@ if __name__ == '__main__':
if ret != 0:
print("[ERROR] simulate returned {0}".format(ret))
- p_x = (c_double * num_atoms)()
- p_y = (c_double * num_atoms)()
- p_z = (c_double * num_atoms)()
- ret = get_atom_positions(handle, p_x, p_y, p_z)
+ qm_p_x = (c_double * num_qm_atoms)()
+ qm_p_y = (c_double * num_qm_atoms)()
+ qm_p_z = (c_double * num_qm_atoms)()
+ mm_p_x = (c_double * num_mm_atoms)()
+ mm_p_y = (c_double * num_mm_atoms)()
+ mm_p_z = (c_double * num_mm_atoms)()
+ ret = get_atom_positions_qmmm(handle, qm_p_x, qm_p_y, qm_p_z,
+ mm_p_x, mm_p_y, mm_p_z)
if ret != 0:
- print("[ERROR] get_atom_positions returned {0}".format(ret))
+ print("[ERROR] get_atom_positions_qmmm returned {0}".format(ret))
- f_x = (c_double * num_atoms)()
- f_y = (c_double * num_atoms)()
- f_z = (c_double * num_atoms)()
- ret = get_atom_forces(handle, f_x, f_y, f_z)
+ qm_f_x = (c_double * num_qm_atoms)()
+ qm_f_y = (c_double * num_qm_atoms)()
+ qm_f_z = (c_double * num_qm_atoms)()
+ mm_f_x = (c_double * num_mm_atoms)()
+ mm_f_y = (c_double * num_mm_atoms)()
+ mm_f_z = (c_double * num_mm_atoms)()
+ ret = get_atom_forces_qmmm(handle, qm_f_x, qm_f_y, qm_f_z, mm_f_x, mm_f_y, mm_f_z)
if ret != 0:
- print("[ERROR] get_atom_forces returned {0}".format(ret))
+ print("[ERROR] get_atom_forces_qmmm returned {0}".format(ret))
- q = (c_double * num_atoms)()
- ret = get_atom_charges(handle, q)
+ qm_q = (c_double * num_qm_atoms)()
+ mm_q = (c_double * num_mm_atoms)()
+ ret = get_atom_charges_qmmm(handle, qm_q, mm_q)
if ret != 0:
- print("[ERROR] get_atom_charges returned {0}".format(ret))
+ print("[ERROR] get_atom_charges_qmmm returned {0}".format(ret))
cleanup(handle)