spuremd.c

    if ( handle != NULL )
    {
        spmd_handle = (spuremd_handle*) handle;

        ret_ = Set_Control_Parameter( keyword, values, spmd_handle->control,
                spmd_handle->out_control );

        if ( ret_ == SUCCESS )
        {
            ret = SPUREMD_SUCCESS;
        }
    }

    return ret;
}


/* Setter for contiguous charge constraints
 *
 * handle: pointer to wrapper struct with top-level data structures
 * num_charge_constraint_contig: num. of contiguous charge constraints for charge model
 * charge_constraint_contig_start: starting atom num. (1-based) of atom group for a charge constraint
 * charge_constraint_contig_end: ending atom num. (1-based) of atom group for a charge constraint
 * charge_constraint_contig_value: charge constraint value for atom group
 *
 * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
 */
int set_contiguous_charge_constraints( const void * const handle,
        int num_charge_constraint_contig, const int * const charge_constraint_contig_start,
        const int * const charge_constraint_contig_end,
        const double * const charge_constraint_contig_value )
{
    int i, ret;
    spuremd_handle *spmd_handle;

    ret = SPUREMD_FAILURE;

    if ( handle != NULL )
    {
        spmd_handle = (spuremd_handle*) handle;

        spmd_handle->system->num_molec_charge_constraints = num_charge_constraint_contig;

        if ( spmd_handle->system->num_molec_charge_constraints
                > spmd_handle->system->max_num_molec_charge_constraints )
        {
            if ( spmd_handle->system->max_num_molec_charge_constraints > 0 )
            {
                sfree( spmd_handle->system->molec_charge_constraints,
                        __FILE__, __LINE__ );
                sfree( spmd_handle->system->molec_charge_constraint_ranges,
                        __FILE__, __LINE__ );
            }

            spmd_handle->system->molec_charge_constraints = smalloc(
                    sizeof(real) * spmd_handle->system->num_molec_charge_constraints,
                    __FILE__, __LINE__ );
            spmd_handle->system->molec_charge_constraint_ranges = smalloc(
                    sizeof(int) * 2 * spmd_handle->system->num_molec_charge_constraints,
                    __FILE__, __LINE__ );

            spmd_handle->system->max_num_molec_charge_constraints
                = spmd_handle->system->num_molec_charge_constraints;
        }

        if ( spmd_handle->system->num_molec_charge_constraints > 0 )
        {
            for ( i = 0; i < spmd_handle->system->num_molec_charge_constraints; ++i )
            {
                spmd_handle->system->molec_charge_constraint_ranges[2 * i] = charge_constraint_contig_start[i];
                spmd_handle->system->molec_charge_constraint_ranges[2 * i + 1] = charge_constraint_contig_end[i];
            }

            for ( i = 0; i < spmd_handle->system->num_molec_charge_constraints; ++i )
            {
                spmd_handle->system->molec_charge_constraints[i] = charge_constraint_contig_value[i];
            }
        }

        ret = SPUREMD_SUCCESS;
    }

    return ret;
}


/* Setter for custom charge constraints
 *
 * handle: pointer to wrapper struct with top-level data structures
 * num_charge_constraint_custom: num. of custom charge constraints for charge model
 * charge_constraint_custom_count: counts for each custom charge constraint
 * charge_constraint_custom_atom_index: atom indices (1-based) for custom charge constraints
 * charge_constraint_custom_coeff: coefficients for custom charge constraints
 * charge_constraint_custom_rhs: right-hand side (RHS) constants for custom charge constraints
 *
 * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
 */
int set_custom_charge_constraints( const void * const handle,
        int num_charge_constraint_custom,
        const int * const charge_constraint_custom_count,
        const int * const charge_constraint_custom_atom_index,
        const double * const charge_constraint_custom_coeff,
        const double * const charge_constraint_custom_rhs )
{
    int i, ret;
    spuremd_handle *spmd_handle;

    ret = SPUREMD_FAILURE;

    if ( handle != NULL )
    {
        spmd_handle = (spuremd_handle*) handle;

        spmd_handle->system->num_custom_charge_constraints = num_charge_constraint_custom;

        if ( spmd_handle->system->num_custom_charge_constraints
                > spmd_handle->system->max_num_custom_charge_constraints )
        {
            if ( spmd_handle->system->max_num_custom_charge_constraints > 0 )
            {
                sfree( spmd_handle->system->custom_charge_constraint_count,
                        __FILE__, __LINE__ );
                sfree( spmd_handle->system->custom_charge_constraint_start,
                        __FILE__, __LINE__ );
                sfree( spmd_handle->system->custom_charge_constraint_rhs,
                        __FILE__, __LINE__ );
            }

            spmd_handle->system->custom_charge_constraint_count = smalloc(
                    sizeof(int) * spmd_handle->system->num_custom_charge_constraints,
                    __FILE__, __LINE__ );
            spmd_handle->system->custom_charge_constraint_start = smalloc(
                    sizeof(int) * (spmd_handle->system->num_custom_charge_constraints + 1),
                    __FILE__, __LINE__ );
            spmd_handle->system->custom_charge_constraint_rhs = smalloc(
                    sizeof(real) * spmd_handle->system->num_custom_charge_constraints,
                    __FILE__, __LINE__ );

            spmd_handle->system->max_num_custom_charge_constraints
                = spmd_handle->system->num_custom_charge_constraints;
        }

        spmd_handle->system->num_custom_charge_constraint_entries = 0;

        if ( spmd_handle->system->num_custom_charge_constraints > 0 )
        {
            for ( i = 0; i < spmd_handle->system->num_custom_charge_constraints; ++i )
            {
                spmd_handle->system->custom_charge_constraint_count[i] = charge_constraint_custom_count[i];
                spmd_handle->system->custom_charge_constraint_start[i] = (i == 0 ? 0 :
                        spmd_handle->system->custom_charge_constraint_start[i - 1] + charge_constraint_custom_count[i - 1]);
                spmd_handle->system->num_custom_charge_constraint_entries += charge_constraint_custom_count[i];
                spmd_handle->system->custom_charge_constraint_rhs[i] = charge_constraint_custom_rhs[i];
            }
            spmd_handle->system->custom_charge_constraint_start[spmd_handle->system->num_custom_charge_constraints]
                = spmd_handle->system->custom_charge_constraint_start[spmd_handle->system->num_custom_charge_constraints - 1]
                + charge_constraint_custom_count[spmd_handle->system->num_custom_charge_constraints - 1];
        }

        if ( spmd_handle->system->num_custom_charge_constraint_entries
                > spmd_handle->system->max_num_custom_charge_constraint_entries )
        {
            if ( spmd_handle->system->max_num_custom_charge_constraint_entries > 0 )
            {
                sfree( spmd_handle->system->custom_charge_constraint_atom_index,
                        __FILE__, __LINE__ );
                sfree( spmd_handle->system->custom_charge_constraint_coeff,
                        __FILE__, __LINE__ );
            }

            spmd_handle->system->custom_charge_constraint_atom_index = smalloc(
                    sizeof(int) * spmd_handle->system->num_custom_charge_constraint_entries,
                    __FILE__, __LINE__ );
            spmd_handle->system->custom_charge_constraint_coeff = smalloc(
                    sizeof(real) * spmd_handle->system->num_custom_charge_constraint_entries,
                    __FILE__, __LINE__ );

            spmd_handle->system->max_num_custom_charge_constraint_entries
                = spmd_handle->system->num_custom_charge_constraint_entries;
        }

        if ( spmd_handle->system->num_custom_charge_constraint_entries > 0 )
        {
            for ( i = 0; i < spmd_handle->system->num_custom_charge_constraint_entries; ++i )
            {
                spmd_handle->system->custom_charge_constraint_atom_index[i] = charge_constraint_custom_atom_index[i];
            }

            for ( i = 0; i < spmd_handle->system->num_custom_charge_constraint_entries; ++i )
            {
                spmd_handle->system->custom_charge_constraint_coeff[i] = charge_constraint_custom_coeff[i];
            }
        }

        ret = SPUREMD_SUCCESS;
    }

    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_symbols: element types for QM atoms
 * qm_pos: coordinates of QM atom positions (consecutively arranged), in Angstroms
 * mm_num_atoms: num. atoms in the MM region
 * mm_symbols: element types for MM atoms
 * mm_pos_q: coordinates and charges of MM atom positions (consecutively arranged), in Angstroms / 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 char * const qm_symbols,
        const double * const qm_pos, int mm_num_atoms, const char * const mm_symbols,
        const double * const mm_pos_q, const double * const sim_box_info,
        const char * const ffield_file, const char * const control_file )
{
    int i;
    char element[3];
    rvec x;
    spuremd_handle *spmd_handle;

    Allocate_Top_Level_Structs( &spmd_handle );
    Initialize_Top_Level_Structs( spmd_handle );

    /* override default */
    spmd_handle->output_enabled = 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, FALSE );

    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, (int) CEIL( SAFE_ZONE * 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 );

    element[2] = '\0';

    for ( i = 0; i < spmd_handle->system->N_qm; ++i )
    {
        x[0] = qm_pos[3 * i];
        x[1] = qm_pos[3 * i + 1];
        x[2] = qm_pos[3 * i + 2];

        Fit_to_Periodic_Box( &spmd_handle->system->box, x );

        spmd_handle->workspace->orig_id[i] = i + 1;
        element[0] = toupper( qm_symbols[2 * i] );
        element[1] = toupper( qm_symbols[2 * i + 1] );
        Trim_Spaces( element, sizeof(element) );
        spmd_handle->system->atoms[i].type = Get_Atom_Type( &spmd_handle->system->reax_param,
                element, sizeof(element) );
        strncpy( spmd_handle->system->atoms[i].name, element,
                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;
            
        /* check for dummy atom */
        if ( strncmp( element, "X\0", 2 ) == 0 )
        {
           spmd_handle->system->atoms[i].is_dummy = TRUE;
        }
        else
        {
            spmd_handle->system->atoms[i].is_dummy = FALSE;            
        }		
    }

    for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
    {
        x[0] = mm_pos_q[4 * (i - spmd_handle->system->N_qm)];
        x[1] = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 1];
        x[2] = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 2];

        Fit_to_Periodic_Box( &spmd_handle->system->box, x );

        spmd_handle->workspace->orig_id[i] = i + 1;
        element[0] = toupper( mm_symbols[2 * (i - spmd_handle->system->N_qm)] );
        element[1] = toupper( mm_symbols[2 * (i - spmd_handle->system->N_qm) + 1] );
        Trim_Spaces( element, sizeof(element) );
        spmd_handle->system->atoms[i].type = Get_Atom_Type( &spmd_handle->system->reax_param,
                element, sizeof(element) );
        strncpy( spmd_handle->system->atoms[i].name, element,
                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_pos_q[4 * (i - spmd_handle->system->N_qm) + 3];
        spmd_handle->system->atoms[i].q_init = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 3];
        spmd_handle->system->atoms[i].qmmm_mask = FALSE;
            
        /* check for dummy atom */
        if ( strncmp( element, "X\0", 2 ) == 0 )
        {
           spmd_handle->system->atoms[i].is_dummy = TRUE;
        }
        else
        {
            spmd_handle->system->atoms[i].is_dummy = FALSE;            
        }		
    }

    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_symbols: element types for QM atoms
 * qm_pos: coordinates of QM atom positions (consecutively arranged), in Angstroms
 * mm_num_atoms: num. atoms in the MM region
 * mm_symbols: element types for MM atoms
 * mm_pos_q: coordinates and charges of MM atom positions (consecutively arranged), in Angstroms / 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 char * const qm_symbols, const double * const qm_pos,
        int mm_num_atoms, const char * const mm_symbols,
        const double * const mm_pos_q, const double * const sim_box_info,
        const char * const ffield_file, const char * const control_file )
{
    int i, ret;
    char element[3];
    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++;

        Read_Input_Files( NULL, ffield_file, control_file,
                spmd_handle->system, spmd_handle->control,
                spmd_handle->data, spmd_handle->workspace,
                spmd_handle->out_control, TRUE );

        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;
        spmd_handle->system->num_molec_charge_constraints = 0;
        spmd_handle->system->num_custom_charge_constraints = 0;
        spmd_handle->system->num_custom_charge_constraint_entries = 0;

        if ( spmd_handle->system->prealloc_allocated == FALSE
                || spmd_handle->system->N > spmd_handle->system->N_max )
        {
            PreAllocate_Space( spmd_handle->system, spmd_handle->control,
                    spmd_handle->workspace, (int) CEIL( SAFE_ZONE * 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 );

        element[2] = '\0';

        for ( i = 0; i < spmd_handle->system->N_qm; ++i )
        {
            x[0] = qm_pos[3 * i];
            x[1] = qm_pos[3 * i + 1];
            x[2] = qm_pos[3 * i + 2];

            Fit_to_Periodic_Box( &spmd_handle->system->box, x );

            spmd_handle->workspace->orig_id[i] = i + 1;
            element[0] = toupper( qm_symbols[2 * i] );
            element[1] = toupper( qm_symbols[2 * i + 1] );
            Trim_Spaces( element, sizeof(element) );
            spmd_handle->system->atoms[i].type = Get_Atom_Type( &spmd_handle->system->reax_param,
                    element, sizeof(element) );
            strncpy( spmd_handle->system->atoms[i].name, element,
                    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;
            
            /* check for dummy atom */
            if ( strncmp( element, "X\0", 2 ) == 0 )
            {
               spmd_handle->system->atoms[i].is_dummy = TRUE;
            }
            else
            {
                spmd_handle->system->atoms[i].is_dummy = FALSE;            
            }		
        }

        for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
        {
            x[0] = mm_pos_q[4 * (i - spmd_handle->system->N_qm)];
            x[1] = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 1];
            x[2] = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 2];

            Fit_to_Periodic_Box( &spmd_handle->system->box, x );

            spmd_handle->workspace->orig_id[i] = i + 1;
            element[0] = toupper( mm_symbols[2 * (i - spmd_handle->system->N_qm)] );
            element[1] = toupper( mm_symbols[2 * (i - spmd_handle->system->N_qm) + 1] );
            Trim_Spaces( element, sizeof(element) );
            spmd_handle->system->atoms[i].type = Get_Atom_Type( &spmd_handle->system->reax_param,
                    element, sizeof(element) );
            strncpy( spmd_handle->system->atoms[i].name, element,
                    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_pos_q[4 * (i - spmd_handle->system->N_qm) + 3];
            spmd_handle->system->atoms[i].q_init = mm_pos_q[4 * (i - spmd_handle->system->N_qm) + 3];
            spmd_handle->system->atoms[i].qmmm_mask = FALSE;
            
            /* check for dummy atom */
            if ( strncmp( element, "X\0", 2 ) == 0 )
            {
               spmd_handle->system->atoms[i].is_dummy = TRUE;
            }
            else
            {
                spmd_handle->system->atoms[i].is_dummy = FALSE;            
            }		
        }

        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: coordinates of QM atom positions (consecutively arranged), in Angstroms (allocated by caller)
 * mm_pos: coordinates of MM atom positions (consecutively arranged), 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,
        double * const mm_pos )
{
    int i, ret;
    spuremd_handle *spmd_handle;

    ret = SPUREMD_FAILURE;

    if ( handle != NULL )
    {
        spmd_handle = (spuremd_handle*) handle;

        if ( qm_pos != NULL )
        {
            for ( i = 0; i < spmd_handle->system->N_qm; ++i )
            {
                qm_pos[3 * i] = spmd_handle->system->atoms[i].x[0];
                qm_pos[3 * i + 1] = spmd_handle->system->atoms[i].x[1];
                qm_pos[3 * i + 2] = spmd_handle->system->atoms[i].x[2];
            }
        }

        if ( mm_pos != NULL )
        {
            for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
            {
                mm_pos[3 * (i - spmd_handle->system->N_qm)] = spmd_handle->system->atoms[i].x[0];
                mm_pos[3 * (i - spmd_handle->system->N_qm) + 1] = spmd_handle->system->atoms[i].x[1];
                mm_pos[3 * (i - spmd_handle->system->N_qm) + 2] = 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: coordinates of QM atom velocities (consecutively arranged), in Angstroms / ps (allocated by caller)
 * mm_vel: coordinates of MM atom velocities (consecutively arranged), 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,
        double * const mm_vel )
{
    int i, ret;
    spuremd_handle *spmd_handle;

    ret = SPUREMD_FAILURE;

    if ( handle != NULL )
    {
        spmd_handle = (spuremd_handle*) handle;

        if ( qm_vel != NULL )
        {
            for ( i = 0; i < spmd_handle->system->N_qm; ++i )
            {
                qm_vel[3 * i] = spmd_handle->system->atoms[i].v[0];
                qm_vel[3 * i + 1] = spmd_handle->system->atoms[i].v[1];
                qm_vel[3 * i + 2] = spmd_handle->system->atoms[i].v[2];
            }
        }

        if ( mm_vel != NULL )
        {
            for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
            {
                mm_vel[3 * (i - spmd_handle->system->N_qm)] = spmd_handle->system->atoms[i].v[0];
                mm_vel[3 * (i - spmd_handle->system->N_qm) + 1] = spmd_handle->system->atoms[i].v[1];
                mm_vel[3 * (i - spmd_handle->system->N_qm) + 2] = 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: coordinates of QM atom forces (consecutively arranged), in Angstroms * Daltons / ps^2 (allocated by caller)
 * mm_f: coordinates of MM atom forces (consecutively arranged), 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,
        double * const mm_f )
{
    int i, ret;
    spuremd_handle *spmd_handle;

    ret = SPUREMD_FAILURE;

    if ( handle != NULL )
    {
        spmd_handle = (spuremd_handle*) handle;

        if ( qm_f != NULL )
        {
            for ( i = 0; i < spmd_handle->system->N_qm; ++i )
            {
                qm_f[3 * i] = spmd_handle->system->atoms[i].f[0];
                qm_f[3 * i + 1] = spmd_handle->system->atoms[i].f[1];
                qm_f[3 * i + 2] = spmd_handle->system->atoms[i].f[2];
            }
        }

        if ( mm_f != NULL )
        {
            for ( i = spmd_handle->system->N_qm; i < spmd_handle->system->N; ++i )
            {
                mm_f[3 * (i - spmd_handle->system->N_qm)] = spmd_handle->system->atoms[i].f[0];
                mm_f[3 * (i - spmd_handle->system->N_qm) + 1] = spmd_handle->system->atoms[i].f[1];
                mm_f[3 * (i - spmd_handle->system->N_qm) + 2] = 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;

        if ( qm_q != NULL )
        {
            for ( i = 0; i < spmd_handle->system->N_qm; ++i )
            {
                qm_q[i] = spmd_handle->system->atoms[i].q;
            }
        }

        if ( mm_q != NULL )
        {
            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