diff --git a/configure.ac b/configure.ac
index 18133ec5aaffafcb1a27c11db895f5535b9e8867..07bcc4539b094206dabc3c013589928615b330d4 100644
--- a/configure.ac
+++ b/configure.ac
@@ -67,6 +67,12 @@ AC_ARG_ENABLE([qmmm],
[enable build for code in QM/MM mode @<:@default: no@:>@])],
[qmmm=${enableval}], [qmmm=no])
+# Build code with Fortran interface in QM/MM mode.
+AC_ARG_ENABLE([qmmm_fortran],
+ [AS_HELP_STRING([--enable-qmmm_fortran],
+ [enable build for code in QM/MM mode @<:@default: no@:>@])],
+ [qmmm_fortran=${enableval}], [qmmm_fortran=no])
+
# Build LAMMPS/reaxc integration code.
AC_ARG_ENABLE([lammps-reaxc],
[AS_HELP_STRING([--enable-lammps-reaxc],
@@ -326,9 +332,13 @@ if test "x${pack_serial_enabled}" = "xyes" || test "x${pack_openmp_enabled}" = "
AC_SUBST([G_LIBS], ["${GSL_LIBS}"])
# Build code in QM/MM mode
- AS_IF([test "x${qmmm}" = "xyes"],
+ AS_IF([test "x${qmmm}" = "xyes" || test "x${qmmm_fortran}" = "xyes"],
[AC_DEFINE([QMMM], [1], [Define to 1 to build PuReMD code in QMMM mode.])])
+ # Build code with Fortran interface in QM/MM mode
+ AS_IF([test "x${qmmm_fortran}" = "xyes"],
+ [AC_DEFINE([QMMM_FORTRAN], [1], [Define to 1 to build PuReMD code in QMMM mode.])])
+
AC_LANG_POP([C])
fi
AM_CONDITIONAL([BUILD_S_OMP], [test "x${pack_serial_enabled}" = "xyes" || test "x${pack_openmp_enabled}" = "xyes"])
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 3d8b8021f864c457df31359abc99610043b2d078..73a87973706eb8c70a7ee20539b677c878a9125a 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -1391,12 +1391,8 @@ static void QEq( reax_system * const system, control_params * const control,
}
#if defined(QMMM)
- for ( int i = 0; i < system->N; ++i )
- {
- workspace->s[0][i] = system->atoms[i].q_init;
- workspace->t[0][i] = system->atoms[i].q_init;
- workspace->mask_qmmm[i] = system->atoms[i].qmmm_mask;
- }
+ fprintf( stderr, "[ERROR] QEq charge method is not supported in QM/MM mode. Use EEM instead. Terminating...\n" );
+ exit( INVALID_INPUT );
#endif
switch ( control->cm_solver_type )
@@ -1412,7 +1408,7 @@ static void QEq( reax_system * const system, control_params * const control,
iters = GMRES_HouseHolder( workspace, control, data, &workspace->H,
workspace->b_s, control->cm_solver_q_err, workspace->s[0], refactor );
iters += GMRES_HouseHolder( workspace, control, data, &workspace->H,
- workspace->b_t, control->cm_solver_q_err, workspace->t[0], 0 );
+ workspace->b_t, control->cm_solver_q_err, workspace->t[0], FALSE );
break;
case CG_S:
@@ -1512,11 +1508,16 @@ static void EE( reax_system * const system, control_params * const control,
}
#if defined(QMMM)
- for ( int i = 0; i < system->N; ++i )
+ for ( int i = 0; i < system->N_qm; ++i )
+ {
+ workspace->mask_qmmm[i] = system->atoms[i].qmmm_mask;
+ }
+ for ( int i = system->N_qm; i < system->N; ++i )
{
workspace->s[0][i] = system->atoms[i].q_init;
workspace->mask_qmmm[i] = system->atoms[i].qmmm_mask;
}
+ workspace->mask_qmmm[system->N_cm - 1] = 1;
#endif
switch ( control->cm_solver_type )
@@ -1631,11 +1632,22 @@ static void ACKS2( reax_system * const system, control_params * const control,
#endif
#if defined(QMMM)
- for ( int i = 0; i < system->N; ++i )
+ /* TODO: further testing needed for QM/MM mode with ACKS2 */
+ for ( int i = 0; i < system->N_qm; ++i )
+ {
+ workspace->mask_qmmm[i] = system->atoms[i].qmmm_mask;
+ }
+ for ( int i = system->N_qm; i < system->N; ++i )
{
workspace->s[0][i] = system->atoms[i].q_init;
workspace->mask_qmmm[i] = system->atoms[i].qmmm_mask;
}
+ for ( int i = system->N; i < 2 * system->N; ++i )
+ {
+ 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;
#endif
switch ( control->cm_solver_type )
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index e850ce3218f6aabfeae08569176ee68a3cd04432..16eb339829d948b619bc6163aae14d7f04d415dc 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -684,7 +684,7 @@ static void Init_Workspace( reax_system *system, control_params *control,
}
#if defined(QMMM)
- workspace->mask_qmmm = smalloc( system->N_max * sizeof( int ),
+ workspace->mask_qmmm = smalloc( system->N_cm_max * sizeof( int ),
"Init_Workspace::workspace->mask_qmmm" );
#endif
diff --git a/sPuReMD/src/nonbonded.c b/sPuReMD/src/nonbonded.c
index 74383c8c033bf05811fdfe30698335ff0b03c43c..807226d532052aea191e3889115d4ea3e6ec9627 100644
--- a/sPuReMD/src/nonbonded.c
+++ b/sPuReMD/src/nonbonded.c
@@ -153,10 +153,6 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
nbr_pj = &far_nbrs->far_nbr_list[pj];
j = nbr_pj->nbr;
-#if defined(QMMM)
- if ( system->atoms[i].qmmm_mask == TRUE || system->atoms[j].qmmm_mask == TRUE )
- {
-#endif
r_ij = nbr_pj->d;
twbp = &system->reax_param.tbp[ system->atoms[i].type ]
[ system->atoms[j].type ];
@@ -243,9 +239,6 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
{
e_core = 0.0;
de_core = 0.0;
- e_vdW = 0.0;
- e_base = 0.0;
- de_base = 0.0;
CEvd = 0.0;
}
#endif
@@ -256,11 +249,11 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
e_base, de_base, e_core, de_core ); fflush( stderr );
#endif
+ /* Coulomb Calculations */
#if defined(QMMM)
- if ( system->atoms[i].qmmm_mask == TRUE && system->atoms[j].qmmm_mask == TRUE )
+ if ( system->atoms[i].qmmm_mask == TRUE || system->atoms[j].qmmm_mask == TRUE )
{
#endif
- /* Coulomb Calculations */
dr3gamij_1 = r_ij * r_ij * r_ij + POW( twbp->gamma, -3.0 );
dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
e_clb = C_ELE * (system->atoms[i].q * system->atoms[j].q) / dr3gamij_3;
@@ -274,8 +267,6 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
}
else
{
- e_clb = 0.0;
- e_ele = 0.0;
de_clb = 0.0;
CEclmb = 0.0;
}
@@ -392,9 +383,6 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
rvec_ScaledAdd( workspace->f_vdw[j], +CEvd, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f_ele[i], -CEclmb, nbr_pj->dvec );
rvec_ScaledAdd( workspace->f_ele[j], +CEclmb, nbr_pj->dvec );
-#endif
-#if defined(QMMM)
- }
#endif
}
}
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index eb70ead6c641e94ef5f0e49f6d6d2640605ff563..dfc00cab8d5358d94c7ad12671934c41640a15ef 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -151,156 +151,6 @@ static void Read_Input_Files( const char * const geo_file,
}
-#if defined(QMMM)
-/* Allocate top-level data structures and parse input files
- * for the first simulation
- *
- * qm_num_atoms: num. atoms in the QM region
- * qm_types: element types for QM atoms
- * qm_pos_x: x-coordinate of QM atom positions, in Angstroms
- * qm_pos_y: y-coordinate of QM atom positions, in Angstroms
- * qm_pos_z: z-coordinate of QM atom positions, in Angstroms
- * mm_num_atoms: num. atoms in the MM region
- * mm_types: element types for MM atoms
- * mm_pos_x: x-coordinate of MM atom positions, in Angstroms
- * mm_pos_y: y-coordinate of MM atom positions, in Angstroms
- * mm_pos_z: z-coordinate of MM atom positions, in Angstroms
- * mm_q: charge of MM atom, in Coulombs
- * sim_box_info: simulation box information, where the entries are
- * - box length per dimension (3 entries)
- * - angles per dimension (3 entries)
- * ffield_file: file containing force field parameters
- * control_file: file containing simulation parameters
- */
-void * setup_qmmm_( int qm_num_atoms, const int * const qm_types,
- const double * const qm_pos_x, const double * const qm_pos_y,
- const double * const qm_pos_z, int mm_num_atoms, const int * const mm_types,
- const double * const mm_pos_x, const double * const mm_pos_y,
- const double * const mm_pos_z, const double * const mm_q,
- const double * const sim_box_info, const char * const ffield_file,
- const char * const control_file )
-{
- int i;
-// char atom_name[9];
- rvec x;
- spuremd_handle *spmd_handle;
-
- /* top-level allocation */
- spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
- "setup::spmd_handle" );
-
- /* second-level allocations */
- spmd_handle->system = smalloc( sizeof(reax_system),
- "Setup::spmd_handle->system" );
- spmd_handle->system->prealloc_allocated = FALSE;
- spmd_handle->system->ffield_params_allocated = FALSE;
- spmd_handle->system->g.allocated = FALSE;
-
- spmd_handle->control = smalloc( sizeof(control_params),
- "Setup::spmd_handle->control" );
-
- spmd_handle->data = smalloc( sizeof(simulation_data),
- "Setup::spmd_handle->data" );
-
- spmd_handle->workspace = smalloc( sizeof(static_storage),
- "Setup::spmd_handle->workspace" );
- spmd_handle->workspace->H.allocated = FALSE;
- spmd_handle->workspace->H_full.allocated = FALSE;
- spmd_handle->workspace->H_sp.allocated = FALSE;
- spmd_handle->workspace->H_p.allocated = FALSE;
- spmd_handle->workspace->H_spar_patt.allocated = FALSE;
- spmd_handle->workspace->H_spar_patt_full.allocated = FALSE;
- spmd_handle->workspace->H_app_inv.allocated = FALSE;
- spmd_handle->workspace->L.allocated = FALSE;
- spmd_handle->workspace->U.allocated = FALSE;
-
- spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
- "Setup::spmd_handle->lists" );
- for ( i = 0; i < LIST_N; ++i )
- {
- spmd_handle->lists[i] = smalloc( sizeof(reax_list),
- "Setup::spmd_handle->lists[i]" );
- spmd_handle->lists[i]->allocated = FALSE;
- }
- spmd_handle->out_control = smalloc( sizeof(output_controls),
- "Setup::spmd_handle->out_control" );
-
- spmd_handle->output_enabled = FALSE;
- spmd_handle->realloc = TRUE;
- 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 = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
-
- PreAllocate_Space( spmd_handle->system, spmd_handle->control,
- spmd_handle->workspace, spmd_handle->system->N );
-
- Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
- sim_box_info[3], sim_box_info[4], sim_box_info[5],
- &spmd_handle->system->box );
-
- for ( i = 0; i < spmd_handle->system->N_qm; ++i )
- {
- x[0] = qm_pos_x[i];
- x[1] = qm_pos_y[i];
- x[2] = qm_pos_z[i];
-
- Fit_to_Periodic_Box( &spmd_handle->system->box, x );
-
- spmd_handle->workspace->orig_id[i] = i + 1;
-// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
-// element, sizeof(element) );
- spmd_handle->system->atoms[i].type = qm_types[i];
-// strncpy( spmd_handle->system->atoms[i].name, atom_name,
-// sizeof(spmd_handle->system->atoms[i].name) - 1 );
-// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
- rvec_Copy( spmd_handle->system->atoms[i].x, x );
- rvec_MakeZero( spmd_handle->system->atoms[i].v );
- rvec_MakeZero( spmd_handle->system->atoms[i].f );
- spmd_handle->system->atoms[i].q = 0.0;
- spmd_handle->system->atoms[i].q_init = 0.0;
-
- spmd_handle->system->atoms[i].qmmm_mask = TRUE;
- }
-
- for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
- {
- x[0] = mm_pos_x[i - spmd_handle->system->N_qm];
- x[1] = mm_pos_y[i - spmd_handle->system->N_qm];
- x[2] = mm_pos_z[i - spmd_handle->system->N_qm];
-
- Fit_to_Periodic_Box( &spmd_handle->system->box, x );
-
- spmd_handle->workspace->orig_id[i] = i + 1;
-// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
-// element, sizeof(element) );
- spmd_handle->system->atoms[i].type = mm_types[i - spmd_handle->system->N_qm];
-// strncpy( spmd_handle->system->atoms[i].name, atom_name,
-// sizeof(spmd_handle->system->atoms[i].name) - 1 );
-// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
- rvec_Copy( spmd_handle->system->atoms[i].x, x );
- rvec_MakeZero( spmd_handle->system->atoms[i].v );
- rvec_MakeZero( spmd_handle->system->atoms[i].f );
- spmd_handle->system->atoms[i].q = mm_q[i - spmd_handle->system->N_qm];
- spmd_handle->system->atoms[i].q_init = mm_q[i - spmd_handle->system->N_qm];
-
- spmd_handle->system->atoms[i].qmmm_mask = FALSE;
- }
-
- Read_Input_Files( NULL, ffield_file, control_file,
- spmd_handle->system, spmd_handle->control,
- spmd_handle->data, spmd_handle->workspace,
- spmd_handle->out_control );
-
- spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
-
- return (void *) spmd_handle;
-}
-#endif
-
-
/* Allocate top-level data structures and parse input files
* for the first simulation
*
@@ -588,42 +438,20 @@ int cleanup( const void * const handle )
}
-#if defined(QMMM)
/* Reset for the next simulation by parsing input files and triggering
* reallocation if more space is needed
*
* handle: pointer to wrapper struct with top-level data structures
- * qm_num_atoms: num. atoms in the QM region
- * qm_types: element types for QM atoms
- * qm_pos_x: x-coordinate of QM atom positions, in Angstroms
- * qm_pos_y: y-coordinate of QM atom positions, in Angstroms
- * qm_pos_z: z-coordinate of QM atom positions, in Angstroms
- * mm_num_atoms: num. atoms in the MM region
- * mm_types: element types for MM atoms
- * mm_pos_x: x-coordinate of MM atom positions, in Angstroms
- * mm_pos_y: y-coordinate of MM atom positions, in Angstroms
- * mm_pos_z: z-coordinate of MM atom positions, in Angstroms
- * mm_q: charge of MM atom, in Coulombs
- * sim_box_info: simulation box information, where the entries are
- * - box length per dimension (3 entries)
- * - angles per dimension (3 entries)
+ * geo_file: file containing geometry info of the structure to simulate
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
-int reset_qmmm_( const void * const handle,
- int qm_num_atoms, const int * const qm_types,
- const double * const qm_pos_x, const double * const qm_pos_y,
- const double * const qm_pos_z,
- int mm_num_atoms, const int * const mm_types,
- const double * const mm_pos_x, const double * const mm_pos_y,
- const double * const mm_pos_z, const double * const mm_q,
- const double * const sim_box_info,
+int reset( const void * const handle, const char * const geo_file,
const char * const ffield_file, const char * const control_file )
{
- int i, ret;
- rvec x;
+ int ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
@@ -642,64 +470,7 @@ 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 = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
-
- PreAllocate_Space( spmd_handle->system, spmd_handle->control,
- spmd_handle->workspace, spmd_handle->system->N );
-
- Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
- sim_box_info[3], sim_box_info[4], sim_box_info[5],
- &spmd_handle->system->box );
-
- for ( i = 0; i < spmd_handle->system->N_qm; ++i )
- {
- x[0] = qm_pos_x[i];
- x[1] = qm_pos_y[i];
- x[2] = qm_pos_z[i];
-
- Fit_to_Periodic_Box( &spmd_handle->system->box, x );
-
- spmd_handle->workspace->orig_id[i] = i + 1;
-// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
-// element, sizeof(element) );
- spmd_handle->system->atoms[i].type = qm_types[i];
-// strncpy( spmd_handle->system->atoms[i].name, atom_name,
-// sizeof(spmd_handle->system->atoms[i].name) - 1 );
-// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
- rvec_Copy( spmd_handle->system->atoms[i].x, x );
- rvec_MakeZero( spmd_handle->system->atoms[i].v );
- rvec_MakeZero( spmd_handle->system->atoms[i].f );
- spmd_handle->system->atoms[i].q = 0.0;
-
- spmd_handle->system->atoms[i].qmmm_mask = TRUE;
- }
-
- for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
- {
- x[0] = mm_pos_x[i - spmd_handle->system->N_qm];
- x[1] = mm_pos_y[i - spmd_handle->system->N_qm];
- x[2] = mm_pos_z[i - spmd_handle->system->N_qm];
-
- Fit_to_Periodic_Box( &spmd_handle->system->box, x );
-
- spmd_handle->workspace->orig_id[i] = i + 1;
-// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
-// element, sizeof(element) );
- spmd_handle->system->atoms[i].type = mm_types[i - spmd_handle->system->N_qm];
-// strncpy( spmd_handle->system->atoms[i].name, atom_name,
-// sizeof(spmd_handle->system->atoms[i].name) - 1 );
-// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
- rvec_Copy( spmd_handle->system->atoms[i].x, x );
- rvec_MakeZero( spmd_handle->system->atoms[i].v );
- rvec_MakeZero( spmd_handle->system->atoms[i].f );
- spmd_handle->system->atoms[i].q = mm_q[i - spmd_handle->system->N_qm];
-
- spmd_handle->system->atoms[i].qmmm_mask = FALSE;
- }
-
- Read_Input_Files( NULL, ffield_file, control_file,
+ Read_Input_Files( geo_file, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control );
@@ -722,23 +493,21 @@ int reset_qmmm_( const void * const handle,
return ret;
}
-#endif
-/* Reset for the next simulation by parsing input files and triggering
- * reallocation if more space is needed
+/* Getter for atom positions
*
* handle: pointer to wrapper struct with top-level data structures
- * geo_file: file containing geometry info of the structure to simulate
- * ffield_file: file containing force field parameters
- * control_file: file containing simulation parameters
+ * pos_x: x-coordinate of atom positions, in Angstroms (allocated by caller)
+ * pos_y: y-coordinate of atom positions, in Angstroms (allocated by caller)
+ * pos_z: z-coordinate of atom positions, in Angstroms (allocated by caller)
*
* returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
-int reset( const void * const handle, const char * const geo_file,
- const char * const ffield_file, const char * const control_file )
+int get_atom_positions( const void * const handle, double * const pos_x,
+ double * const pos_y, double * const pos_z )
{
- int ret;
+ int i, ret;
spuremd_handle *spmd_handle;
ret = SPUREMD_FAILURE;
@@ -747,240 +516,7 @@ int reset( const void * const handle, const char * const geo_file,
{
spmd_handle = (spuremd_handle*) handle;
- /* close files used in previous simulation */
- if ( spmd_handle->output_enabled == TRUE )
- {
- Finalize_Out_Controls( spmd_handle->system, spmd_handle->control,
- spmd_handle->workspace, spmd_handle->out_control );
- }
-
- spmd_handle->realloc = FALSE;
- spmd_handle->data->sim_id++;
-
- Read_Input_Files( geo_file, ffield_file, control_file,
- spmd_handle->system, spmd_handle->control,
- spmd_handle->data, spmd_handle->workspace,
- spmd_handle->out_control );
-
- if ( spmd_handle->system->N > spmd_handle->system->N_max )
- {
- /* deallocate everything which needs more space
- * (i.e., structures whose space is a function of the number of atoms),
- * except for data structures allocated while parsing input files */
- Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
- spmd_handle->output_enabled, TRUE );
-
- spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
- spmd_handle->realloc = TRUE;
- }
-
- ret = SPUREMD_SUCCESS;
- }
-
- return ret;
-}
-
-
-#if defined(QMMM)
-/* 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;
-}
-#endif
-
-
-/* Getter for atom positions
- *
- * handle: pointer to wrapper struct with top-level data structures
- * pos_x: x-coordinate of atom positions, in Angstroms (allocated by caller)
- * pos_y: y-coordinate of atom positions, in Angstroms (allocated by caller)
- * pos_z: z-coordinate of atom positions, in Angstroms (allocated by caller)
- *
- * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
- */
-int get_atom_positions( const void * const handle, double * const pos_x,
- double * const pos_y, double * const 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; ++i )
+ for ( i = 0; i < spmd_handle->system->N; ++i )
{
pos_x[i] = spmd_handle->system->atoms[i].x[0];
pos_y[i] = spmd_handle->system->atoms[i].x[1];
@@ -1188,3 +724,751 @@ int set_control_parameter( const void * const handle, const char * const keyword
return ret;
}
+
+
+#if defined(QMMM)
+/* Allocate top-level data structures and parse input files
+ * for the first simulation
+ *
+ * qm_num_atoms: num. atoms in the QM region
+ * qm_types: element types for QM atoms
+ * qm_pos_x: x-coordinate of QM atom positions, in Angstroms
+ * qm_pos_y: y-coordinate of QM atom positions, in Angstroms
+ * qm_pos_z: z-coordinate of QM atom positions, in Angstroms
+ * mm_num_atoms: num. atoms in the MM region
+ * mm_types: element types for MM atoms
+ * mm_pos_x: x-coordinate of MM atom positions, in Angstroms
+ * mm_pos_y: y-coordinate of MM atom positions, in Angstroms
+ * mm_pos_z: z-coordinate of MM atom positions, in Angstroms
+ * mm_q: charge of MM atom, in Coulombs
+ * sim_box_info: simulation box information, where the entries are
+ * - box length per dimension (3 entries)
+ * - angles per dimension (3 entries)
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ */
+void * setup_qmmm( int qm_num_atoms, const int * const qm_types,
+ const double * const qm_pos_x, const double * const qm_pos_y,
+ const double * const qm_pos_z, int mm_num_atoms, const int * const mm_types,
+ const double * const mm_pos_x, const double * const mm_pos_y,
+ const double * const mm_pos_z, const double * const mm_q,
+ const double * const sim_box_info, const char * const ffield_file,
+ const char * const control_file )
+{
+ int i;
+// char atom_name[9];
+ rvec x;
+ spuremd_handle *spmd_handle;
+
+ /* top-level allocation */
+ spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
+ "setup::spmd_handle" );
+
+ /* second-level allocations */
+ spmd_handle->system = smalloc( sizeof(reax_system),
+ "Setup::spmd_handle->system" );
+ spmd_handle->system->prealloc_allocated = FALSE;
+ spmd_handle->system->ffield_params_allocated = FALSE;
+ spmd_handle->system->g.allocated = FALSE;
+
+ spmd_handle->control = smalloc( sizeof(control_params),
+ "Setup::spmd_handle->control" );
+
+ spmd_handle->data = smalloc( sizeof(simulation_data),
+ "Setup::spmd_handle->data" );
+
+ spmd_handle->workspace = smalloc( sizeof(static_storage),
+ "Setup::spmd_handle->workspace" );
+ spmd_handle->workspace->H.allocated = FALSE;
+ spmd_handle->workspace->H_full.allocated = FALSE;
+ spmd_handle->workspace->H_sp.allocated = FALSE;
+ spmd_handle->workspace->H_p.allocated = FALSE;
+ spmd_handle->workspace->H_spar_patt.allocated = FALSE;
+ spmd_handle->workspace->H_spar_patt_full.allocated = FALSE;
+ spmd_handle->workspace->H_app_inv.allocated = FALSE;
+ spmd_handle->workspace->L.allocated = FALSE;
+ spmd_handle->workspace->U.allocated = FALSE;
+
+ spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
+ "Setup::spmd_handle->lists" );
+ for ( i = 0; i < LIST_N; ++i )
+ {
+ spmd_handle->lists[i] = smalloc( sizeof(reax_list),
+ "Setup::spmd_handle->lists[i]" );
+ spmd_handle->lists[i]->allocated = FALSE;
+ }
+ spmd_handle->out_control = smalloc( sizeof(output_controls),
+ "Setup::spmd_handle->out_control" );
+
+ spmd_handle->output_enabled = FALSE;
+ spmd_handle->realloc = TRUE;
+ 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 = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
+
+ PreAllocate_Space( spmd_handle->system, spmd_handle->control,
+ spmd_handle->workspace, spmd_handle->system->N );
+
+ Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
+ sim_box_info[3], sim_box_info[4], sim_box_info[5],
+ &spmd_handle->system->box );
+
+ for ( i = 0; i < spmd_handle->system->N_qm; ++i )
+ {
+ x[0] = qm_pos_x[i];
+ x[1] = qm_pos_y[i];
+ x[2] = qm_pos_z[i];
+
+ Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+ spmd_handle->workspace->orig_id[i] = i + 1;
+// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+// element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = qm_types[i];
+// strncpy( spmd_handle->system->atoms[i].name, atom_name,
+// sizeof(spmd_handle->system->atoms[i].name) - 1 );
+// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ rvec_Copy( spmd_handle->system->atoms[i].x, x );
+ rvec_MakeZero( spmd_handle->system->atoms[i].v );
+ rvec_MakeZero( spmd_handle->system->atoms[i].f );
+ spmd_handle->system->atoms[i].q = 0.0;
+ spmd_handle->system->atoms[i].q_init = 0.0;
+
+ spmd_handle->system->atoms[i].qmmm_mask = TRUE;
+ }
+
+ for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
+ {
+ x[0] = mm_pos_x[i - spmd_handle->system->N_qm];
+ x[1] = mm_pos_y[i - spmd_handle->system->N_qm];
+ x[2] = mm_pos_z[i - spmd_handle->system->N_qm];
+
+ Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+ spmd_handle->workspace->orig_id[i] = i + 1;
+// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+// element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = mm_types[i - spmd_handle->system->N_qm];
+// strncpy( spmd_handle->system->atoms[i].name, atom_name,
+// sizeof(spmd_handle->system->atoms[i].name) - 1 );
+// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ rvec_Copy( spmd_handle->system->atoms[i].x, x );
+ rvec_MakeZero( spmd_handle->system->atoms[i].v );
+ rvec_MakeZero( spmd_handle->system->atoms[i].f );
+ spmd_handle->system->atoms[i].q = mm_q[i - spmd_handle->system->N_qm];
+ spmd_handle->system->atoms[i].q_init = mm_q[i - spmd_handle->system->N_qm];
+
+ spmd_handle->system->atoms[i].qmmm_mask = FALSE;
+ }
+
+ Read_Input_Files( NULL, ffield_file, control_file,
+ spmd_handle->system, spmd_handle->control,
+ spmd_handle->data, spmd_handle->workspace,
+ spmd_handle->out_control );
+
+ spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
+
+ return (void *) spmd_handle;
+}
+
+
+/* Reset for the next simulation by parsing input files and triggering
+ * reallocation if more space is needed
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * qm_num_atoms: num. atoms in the QM region
+ * qm_types: element types for QM atoms
+ * qm_pos_x: x-coordinate of QM atom positions, in Angstroms
+ * qm_pos_y: y-coordinate of QM atom positions, in Angstroms
+ * qm_pos_z: z-coordinate of QM atom positions, in Angstroms
+ * mm_num_atoms: num. atoms in the MM region
+ * mm_types: element types for MM atoms
+ * mm_pos_x: x-coordinate of MM atom positions, in Angstroms
+ * mm_pos_y: y-coordinate of MM atom positions, in Angstroms
+ * mm_pos_z: z-coordinate of MM atom positions, in Angstroms
+ * mm_q: charge of MM atom, in Coulombs
+ * sim_box_info: simulation box information, where the entries are
+ * - box length per dimension (3 entries)
+ * - angles per dimension (3 entries)
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int reset_qmmm( const void * const handle,
+ int qm_num_atoms, const int * const qm_types,
+ const double * const qm_pos_x, const double * const qm_pos_y,
+ const double * const qm_pos_z,
+ int mm_num_atoms, const int * const mm_types,
+ const double * const mm_pos_x, const double * const mm_pos_y,
+ const double * const mm_pos_z, const double * const mm_q,
+ const double * const sim_box_info,
+ const char * const ffield_file, const char * const control_file )
+{
+ int i, ret;
+ rvec x;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ /* close files used in previous simulation */
+ if ( spmd_handle->output_enabled == TRUE )
+ {
+ Finalize_Out_Controls( spmd_handle->system, spmd_handle->control,
+ spmd_handle->workspace, spmd_handle->out_control );
+ }
+
+ 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 = spmd_handle->system->N_qm + spmd_handle->system->N_mm;
+
+ PreAllocate_Space( spmd_handle->system, spmd_handle->control,
+ spmd_handle->workspace, spmd_handle->system->N );
+
+ Setup_Box( sim_box_info[0], sim_box_info[1], sim_box_info[2],
+ sim_box_info[3], sim_box_info[4], sim_box_info[5],
+ &spmd_handle->system->box );
+
+ for ( i = 0; i < spmd_handle->system->N_qm; ++i )
+ {
+ x[0] = qm_pos_x[i];
+ x[1] = qm_pos_y[i];
+ x[2] = qm_pos_z[i];
+
+ Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+ spmd_handle->workspace->orig_id[i] = i + 1;
+// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+// element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = qm_types[i];
+// strncpy( spmd_handle->system->atoms[i].name, atom_name,
+// sizeof(spmd_handle->system->atoms[i].name) - 1 );
+// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ rvec_Copy( spmd_handle->system->atoms[i].x, x );
+ rvec_MakeZero( spmd_handle->system->atoms[i].v );
+ rvec_MakeZero( spmd_handle->system->atoms[i].f );
+ spmd_handle->system->atoms[i].q = 0.0;
+
+ spmd_handle->system->atoms[i].qmmm_mask = TRUE;
+ }
+
+ for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
+ {
+ x[0] = mm_pos_x[i - spmd_handle->system->N_qm];
+ x[1] = mm_pos_y[i - spmd_handle->system->N_qm];
+ x[2] = mm_pos_z[i - spmd_handle->system->N_qm];
+
+ Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+ spmd_handle->workspace->orig_id[i] = i + 1;
+// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+// element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = mm_types[i - spmd_handle->system->N_qm];
+// strncpy( spmd_handle->system->atoms[i].name, atom_name,
+// sizeof(spmd_handle->system->atoms[i].name) - 1 );
+// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ rvec_Copy( spmd_handle->system->atoms[i].x, x );
+ rvec_MakeZero( spmd_handle->system->atoms[i].v );
+ rvec_MakeZero( spmd_handle->system->atoms[i].f );
+ spmd_handle->system->atoms[i].q = mm_q[i - spmd_handle->system->N_qm];
+
+ spmd_handle->system->atoms[i].qmmm_mask = FALSE;
+ }
+
+ Read_Input_Files( NULL, ffield_file, control_file,
+ spmd_handle->system, spmd_handle->control,
+ spmd_handle->data, spmd_handle->workspace,
+ spmd_handle->out_control );
+
+ if ( spmd_handle->system->N > spmd_handle->system->N_max )
+ {
+ /* deallocate everything which needs more space
+ * (i.e., structures whose space is a function of the number of atoms),
+ * except for data structures allocated while parsing input files */
+ Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
+ spmd_handle->output_enabled, TRUE );
+
+ spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
+ spmd_handle->realloc = TRUE;
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
+/* 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;
+}
+#endif
+
+
+#if defined(QMMM_FORTRAN)
+/* Allocate top-level data structures and parse input files
+ * for the first simulation
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * qm_num_atoms: num. atoms in the QM region
+ * qm_types: element types for QM atoms
+ * qm_pos_x: x-coordinate of QM atom positions, in Angstroms
+ * qm_pos_y: y-coordinate of QM atom positions, in Angstroms
+ * qm_pos_z: z-coordinate of QM atom positions, in Angstroms
+ * mm_num_atoms: num. atoms in the MM region
+ * mm_types: element types for MM atoms
+ * mm_pos_x: x-coordinate of MM atom positions, in Angstroms
+ * mm_pos_y: y-coordinate of MM atom positions, in Angstroms
+ * mm_pos_z: z-coordinate of MM atom positions, in Angstroms
+ * mm_q: charge of MM atom, in Coulombs
+ * sim_box_info: simulation box information, where the entries are
+ * - box length per dimension (3 entries)
+ * - angles per dimension (3 entries)
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ */
+void setup_qmmm_( void * handle, const int * const qm_num_atoms, const int * const qm_types,
+ const double * const qm_pos_x, const double * const qm_pos_y,
+ const double * const qm_pos_z, const int * const mm_num_atoms, const int * const mm_types,
+ const double * const mm_pos_x, const double * const mm_pos_y,
+ const double * const mm_pos_z, const double * const mm_q,
+ const double * const sim_box_info, const char * const ffield_file,
+ const char * const control_file )
+{
+ handle = setup_qmmm( *qm_num_atoms, qm_types, qm_pos_x, qm_pos_y, qm_pos_z,
+ *mm_num_atoms, mm_types, mm_pos_x, mm_pos_y, mm_pos_z, mm_q,
+ sim_box_info, ffield_file, control_file );
+}
+
+
+/* Reset for the next simulation by parsing input files and triggering
+ * reallocation if more space is needed
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * qm_num_atoms: num. atoms in the QM region
+ * qm_types: element types for QM atoms
+ * qm_pos_x: x-coordinate of QM atom positions, in Angstroms
+ * qm_pos_y: y-coordinate of QM atom positions, in Angstroms
+ * qm_pos_z: z-coordinate of QM atom positions, in Angstroms
+ * mm_num_atoms: num. atoms in the MM region
+ * mm_types: element types for MM atoms
+ * mm_pos_x: x-coordinate of MM atom positions, in Angstroms
+ * mm_pos_y: y-coordinate of MM atom positions, in Angstroms
+ * mm_pos_z: z-coordinate of MM atom positions, in Angstroms
+ * mm_q: charge of MM atom, in Coulombs
+ * sim_box_info: simulation box information, where the entries are
+ * - box length per dimension (3 entries)
+ * - angles per dimension (3 entries)
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ */
+void reset_qmmm_( const void * const handle,
+ const int * const qm_num_atoms, const int * const qm_types,
+ const double * const qm_pos_x, const double * const qm_pos_y,
+ const double * const qm_pos_z,
+ const int * const mm_num_atoms, const int * const mm_types,
+ const double * const mm_pos_x, const double * const mm_pos_y,
+ const double * const mm_pos_z, const double * const mm_q,
+ const double * const sim_box_info,
+ const char * const ffield_file, const char * const control_file )
+{
+ int ret;
+
+ ret = reset_qmmm( handle, *qm_num_atoms, qm_types, qm_pos_x, qm_pos_y, qm_pos_z,
+ *mm_num_atoms, mm_types, mm_pos_x, mm_pos_y, mm_pos_z, mm_q,
+ sim_box_info, ffield_file, control_file );
+
+ if ( ret != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+
+
+/* Run the simulation according to the prescribed parameters
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ */
+void simulate_( const void * const handle )
+{
+ int ret;
+
+ ret = simulate( handle );
+
+ if ( ret != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+
+
+/* Deallocate all data structures post-simulation
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ */
+void cleanup_( const void * const handle )
+{
+ int ret;
+
+ ret = cleanup( handle );
+
+ if ( ret != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+
+
+/* Setter for writing output to files
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * enabled: TRUE enables writing output to files, FALSE otherwise
+ */
+void set_output_enabled_( const void * const handle, const int enabled )
+{
+ int ret;
+
+ ret = set_output_enabled( handle, enabled );
+
+ if ( ret != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+
+
+/* Setter for simulation parameter values as defined in the input control file
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * control_keyword: keyword from the control file to set the value for
+ * control_value: value to set
+ */
+void set_control_parameter_( const void * const handle, const char * const keyword,
+ const char ** const values )
+{
+ int ret;
+
+ ret = set_control_parameter( handle, keyword, values );
+
+ if ( ret != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+
+
+/* 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)
+ */
+void 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 ret;
+
+ ret = get_atom_positions_qmmm( handle, qm_pos_x, qm_pos_y, qm_pos_z,
+ mm_pos_x, mm_pos_y, mm_pos_z );
+
+ if ( ret != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+
+
+/* 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)
+ */
+void 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 ret;
+
+ ret = get_atom_velocities_qmmm( handle, qm_vel_x, qm_vel_y, qm_vel_z,
+ mm_vel_x, mm_vel_y, mm_vel_z );
+
+ if ( ret != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+
+
+/* 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)
+ */
+void 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 ret;
+
+ 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 != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+
+
+/* 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
+ */
+void get_atom_charges_qmmm_( const void * const handle, double * const qm_q,
+ double * const mm_q )
+{
+ int ret;
+
+ ret = get_atom_charges_qmmm( handle, qm_q, mm_q );
+
+ if ( ret != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+
+
+/* Getter for system energies
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * e_pot: system potential energy, in kcal / mol (reference from caller)
+ * e_kin: system kinetic energy, in kcal / mol (reference from caller)
+ * e_tot: system total energy, in kcal / mol (reference from caller)
+ * t_scalar: temperature scalar, in K (reference from caller)
+ * vol: volume of the simulation box, in Angstroms^3 (reference from caller)
+ * pres: average pressure, in K (reference from caller)
+ */
+void get_system_info_( const void * const handle, double * const e_pot,
+ double * const e_kin, double * const e_tot, double * const temp,
+ double * const vol, double * const pres )
+{
+ int ret;
+
+ ret = get_system_info( handle, e_pot, e_kin, e_tot, temp, vol, pres );
+
+ if ( ret != SPUREMD_SUCCESS )
+ {
+ /* TODO: pass errors via another mechanism */
+ ;
+ }
+}
+#endif
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index ea021f8a771bec2ef6f21bcfafb49aca17ae4535..30251c44f4a2eaf60b111514312efdcf8318d682 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -33,16 +33,6 @@
extern "C" {
#endif
-#if defined(QMMM)
-void * setup_qmmm_( int, const int * const,
- const double * const, const double * const,
- const double * const, int, const int * const,
- const double * const, const double * const,
- const double * const, const double * const,
- const double * const,
- const char * const, const char * const );
-#endif
-
void * setup( const char * const, const char * const,
const char * const );
@@ -52,55 +42,106 @@ int simulate( const void * const );
int cleanup( const void * const );
+int reset( const void * const, const char * const,
+ const char * const, const char * const );
+
+int get_atom_positions( const void * const, double * const,
+ double * const, double * const );
+
+int get_atom_velocities( const void * const, double * const,
+ double * const, double * const );
+
+int get_atom_forces( const void * const, double * const,
+ double * const, double * const );
+
+int get_atom_charges( const void * const, double * const );
+
+int get_system_info( const void * const, double * const,
+ double * const, double * const, double * const,
+ double * const, double * const );
+
+int set_output_enabled( const void * const, const int );
+
+int set_control_parameter( const void * const, const char * const,
+ const char ** const );
+
#if defined(QMMM)
-int reset_qmmm_( const void * const, int,
- const int * const,
+void * setup_qmmm( int, const int * const,
const double * const, const double * const,
const double * const, int, const int * const,
const double * const, const double * const,
const double * const, const double * const,
const double * const,
const char * const, const char * const );
-#endif
-int reset( const void * const, const char * const,
+int reset_qmmm( const void * const, int,
+ const int * const,
+ const double * const, const double * const,
+ const double * const, int, const int * const,
+ const double * const, const double * const,
+ const double * const, const double * const,
+ const double * const,
const char * const, const char * const );
-#if defined(QMMM)
-int get_atom_positions_qmmm_( const void * const, double * 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,
+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,
+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_charges_qmmm( const void * const, double * const, double * const );
#endif
-int get_atom_positions( const void * const, double * const,
- double * const, double * const );
+#if defined(QMMM_FORTRAN)
+void setup_qmmm_( void *, const int * const, const int * const,
+ const double * const, const double * const,
+ const double * const, const int * const, const int * const,
+ const double * const, const double * const,
+ const double * const, const double * const,
+ const double * const,
+ const char * const, const char * const );
-int get_atom_velocities( const void * const, double * const,
- double * const, double * const );
+void reset_qmmm_( const void * const, const int * const, const int * const,
+ const double * const, const double * const,
+ const double * const, const int * const, const int * const,
+ const double * const, const double * const,
+ const double * const, const double * const,
+ const double * const,
+ const char * const, const char * const );
-int get_atom_forces( const void * const, double * const,
+void simulate_( const void * const );
+
+void cleanup_( const void * const );
+
+void set_output_enabled_( const void * const, const int );
+
+void set_control_parameter_( const void * const, const char * const,
+ const char ** const );
+
+void get_atom_positions_qmmm_( const void * const, double * const,
+ double * const, double * const, double * const,
double * const, double * const );
-int get_atom_charges( const void * const, double * const );
+void get_atom_velocities_qmmm_( const void * const, double * const,
+ double * const, double * const, double * const,
+ double * const, double * const );
-int get_system_info( const void * const, double * const,
+void get_atom_forces_qmmm_( const void * const, double * const,
double * const, double * const, double * const,
double * const, double * const );
-int set_output_enabled( const void * const, const int );
+void get_atom_charges_qmmm_( const void * const, double * const, double * const );
-int set_control_parameter( const void * const, const char * const,
- const char ** const );
+void get_system_info_( const void * const, double * const,
+ double * const, double * const, double * const,
+ double * const, double * const );
+#endif
#if defined(__cplusplus)
}