spuremd.c

#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)
 * num_charge_constraints: num. of charge constraints for charge model
 * charge_constraint_start: starting atom num. (1-based) of atom group for a charge constraint
 * charge_constraint_end: ending atom num. (1-based) of atom group for a charge constraint
 * charge_constraint_value: charge constraint value for atom group
 * 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,
        int num_charge_constraints, const int * const charge_constraint_start,
        const int * const charge_constraint_end, const double * const charge_constraint_value,
        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->max_num_molec_charge_constraints = num_charge_constraints;
    spmd_handle->system->num_molec_charge_constraints = num_charge_constraints;

    if ( spmd_handle->system->num_molec_charge_constraints > 0 )
    {
        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__ );

        for ( i = 0; i < spmd_handle->system->num_molec_charge_constraints; ++i )
        {
            spmd_handle->system->molec_charge_constraint_ranges[2 * i] = charge_constraint_start[i];
            spmd_handle->system->molec_charge_constraint_ranges[2 * i + 1] = charge_constraint_end[i];
            spmd_handle->system->molec_charge_constraints[i] = charge_constraint_value[i];
        }
    }

    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;
    /* note: assign here to avoid compiler warning
     * of uninitialized usage in PreAllocate_Space */
    spmd_handle->system->N_max = 0;

    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)
 * num_charge_constraints: num. of charge constraints for charge model
 * charge_constraint_start: starting atom num. (1-based) of atom group for a charge constraint
 * charge_constraint_end: ending atom num. (1-based) of atom group for a charge constraint
 * charge_constraint_value: charge constraint value for atom group
 * 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,
        int num_charge_constraints, const int * const charge_constraint_start,
        const int * const charge_constraint_end, const double * const charge_constraint_value,
        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->num_molec_charge_constraints = num_charge_constraints;

        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_start[i];
                spmd_handle->system->molec_charge_constraint_ranges[2 * i + 1] = charge_constraint_end[i];
                spmd_handle->system->molec_charge_constraints[i] = charge_constraint_value[i];
            }
        }

        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;

        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