init_md.c

    strcpy( temp, control->sim_name );
    strcat( temp, ".f3body" );
    out_control->f3body = fopen( temp, "w" );

    /* open hydrogen bond forces file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".fhb" );
    out_control->fhb = fopen( temp, "w" );

    /* open torsion forces file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".f4body" );
    out_control->f4body = fopen( temp, "w" );

    /* open nonbonded forces file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".fnonb" );
    out_control->fnonb = fopen( temp, "w" );

    /* open total force file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".ftot" );
    out_control->ftot = fopen( temp, "w" );

    /* open coulomb forces file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".ftot2" );
    out_control->ftot2 = fopen( temp, "w" );
#endif


    /* Error handling */
    /* if ( out_control->out == NULL || out_control->pot == NULL ||
       out_control->log == NULL || out_control->mol == NULL ||
       out_control->dpl == NULL || out_control->drft == NULL ||
       out_control->pdb == NULL )
       {
       fprintf( stderr, "FILE OPEN ERROR. TERMINATING..." );
       exit( CANNOT_OPEN_FILE );
       }*/

#undef TEMP_SIZE
}


void Initialize( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace, reax_list **lists,
        output_controls *out_control, evolve_function *Evolve,
        interaction_function *interaction_functions, const int output_enabled )
{
#if defined(DEBUG)
    real start, end;
#endif

#ifdef _OPENMP
    #pragma omp parallel default(none) shared(control)
    {
        #pragma omp single
        control->num_threads = omp_get_num_threads( );
    }
#else
    control->num_threads = 1;
#endif

    Randomize( );

    Init_System( system, control, data );

    Init_Simulation_Data( system, control, data, out_control, Evolve );

    Init_Workspace( system, control, workspace );

    Init_Lists( system, control, data, workspace, lists, out_control );

    if ( output_enabled == TRUE )
    {
        Init_Out_Controls( system, control, workspace, out_control );
    }

    /* These are done in forces.c, only forces.c can see all those functions */
    Init_Bonded_Force_Functions( control, interaction_functions );

#ifdef TEST_FORCES
    Init_Force_Test_Functions( );
#endif

    if ( control->tabulate )
    {
#if defined(DEBUG)
        start = Get_Time( );
#endif

        Make_LR_Lookup_Table( system, control, workspace );

#if defined(DEBUG)
        end = Get_Timing_Info( start );

        fprintf( stderr, "Time for LR Lookup Table calculation is %f \n", end );
#endif
    }

#if defined(DEBUG_FOCUS)
    fprintf( stderr, "data structures have been initialized...\n" );
#endif
}


void Finalize_System( reax_system *system, control_params *control,
        simulation_data *data )
{
    int i, j, k;
    reax_interaction *reax;

    reax = &( system->reaxprm );

    Finalize_Grid( system );

    sfree( reax->gp.l, "Finalize_System::reax->gp.l" );

    for ( i = 0; i < reax->num_atom_types; i++ )
    {
        for ( j = 0; j < reax->num_atom_types; j++ )
        {
            for ( k = 0; k < reax->num_atom_types; k++ )
            {
                sfree( reax->fbp[i][j][k], "Finalize_System::reax->fbp[i][j][k]" );
            }

            sfree( reax->thbp[i][j], "Finalize_System::reax->thbp[i][j]" );
            sfree( reax->hbp[i][j], "Finalize_System::reax->hbp[i][j]" );
            sfree( reax->fbp[i][j], "Finalize_System::reax->fbp[i][j]" );
        }

        sfree( reax->tbp[i], "Finalize_System::reax->tbp[i]" );
        sfree( reax->thbp[i], "Finalize_System::reax->thbp[i]" );
        sfree( reax->hbp[i], "Finalize_System::reax->hbp[i]" );
        sfree( reax->fbp[i], "Finalize_System::reax->fbp[i]" );
    }

    sfree( reax->sbp, "Finalize_System::reax->sbp" );
    sfree( reax->tbp, "Finalize_System::reax->tbp" );
    sfree( reax->thbp, "Finalize_System::reax->thbp" );
    sfree( reax->hbp, "Finalize_System::reax->hbp" );
    sfree( reax->fbp, "Finalize_System::reax->fbp" );

    sfree( system->atoms, "Finalize_System::system->atoms" );
}


void Finalize_Simulation_Data( reax_system *system, control_params *control,
        simulation_data *data, output_controls *out_control )
{
}


void Finalize_Workspace( reax_system *system, control_params *control,
        static_storage *workspace )
{
    int i;

    sfree( workspace->hbond_index, "Finalize_Workspace::workspace->hbond_index" );
    sfree( workspace->total_bond_order, "Finalize_Workspace::workspace->total_bond_order" );
    sfree( workspace->Deltap, "Finalize_Workspace::workspace->Deltap" );
    sfree( workspace->Deltap_boc, "Finalize_Workspace::workspace->Deltap_boc" );
    sfree( workspace->dDeltap_self, "Finalize_Workspace::workspace->dDeltap_self" );
    sfree( workspace->Delta, "Finalize_Workspace::workspace->Delta" );
    sfree( workspace->Delta_lp, "Finalize_Workspace::workspace->Delta_lp" );
    sfree( workspace->Delta_lp_temp, "Finalize_Workspace::workspace->Delta_lp_temp" );
    sfree( workspace->dDelta_lp, "Finalize_Workspace::workspace->dDelta_lp" );
    sfree( workspace->dDelta_lp_temp, "Finalize_Workspace::workspace->dDelta_lp_temp" );
    sfree( workspace->Delta_e, "Finalize_Workspace::workspace->Delta_e" );
    sfree( workspace->Delta_boc, "Finalize_Workspace::workspace->Delta_boc" );
    sfree( workspace->nlp, "Finalize_Workspace::workspace->nlp" );
    sfree( workspace->nlp_temp, "Finalize_Workspace::workspace->nlp_temp" );
    sfree( workspace->Clp, "Finalize_Workspace::workspace->Clp" );
    sfree( workspace->CdDelta, "Finalize_Workspace::workspace->CdDelta" );
    sfree( workspace->vlpex, "Finalize_Workspace::workspace->vlpex" );

    Deallocate_Matrix( workspace->H );
    Deallocate_Matrix( workspace->H_sp );
    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
            control->cm_solver_pre_comp_type == ILU_PAR_PC ||
            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
    {
        Deallocate_Matrix( workspace->L );
        Deallocate_Matrix( workspace->U );
    }
    if ( control->cm_solver_pre_comp_type == SAI_PC )
    {
        Deallocate_Matrix( workspace->H_full );
        Deallocate_Matrix( workspace->H_spar_patt );
        Deallocate_Matrix( workspace->H_spar_patt_full );
        Deallocate_Matrix( workspace->H_app_inv );
    }
    if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
    {
        Deallocate_Matrix( workspace->H_full );
        Deallocate_Matrix( workspace->H_p );
    }

    for ( i = 0; i < 5; ++i )
    {
        sfree( workspace->s[i], "Finalize_Workspace::workspace->s[i]" );
        sfree( workspace->t[i], "Finalize_Workspace::workspace->t[i]" );
    }

    if ( control->cm_solver_pre_comp_type == DIAG_PC )
    {
        sfree( workspace->Hdia_inv, "Finalize_Workspace::workspace->Hdia_inv" );
    }
    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
    {
        sfree( workspace->droptol, "Finalize_Workspace::workspace->droptol" );
    }
    //TODO: check if unused
    //sfree( workspace->w, "Finalize_Workspace::workspace->w" );
    //TODO: check if unused
    sfree( workspace->b, "Finalize_Workspace::workspace->b" );
    sfree( workspace->b_s, "Finalize_Workspace::workspace->b_s" );
    sfree( workspace->b_t, "Finalize_Workspace::workspace->b_t" );
    sfree( workspace->b_prc, "Finalize_Workspace::workspace->b_prc" );
    sfree( workspace->b_prm, "Finalize_Workspace::workspace->b_prm" );
    sfree( workspace->s, "Finalize_Workspace::workspace->s" );
    sfree( workspace->t, "Finalize_Workspace::workspace->t" );

    switch ( control->cm_solver_type )
    {
        /* GMRES storage */
        case GMRES_S:
        case GMRES_H_S:
            for ( i = 0; i < control->cm_solver_restart + 1; ++i )
            {
                sfree( workspace->h[i], "Finalize_Workspace::workspace->h[i]" );
                sfree( workspace->rn[i], "Finalize_Workspace::workspace->rn[i]" );
                sfree( workspace->v[i], "Finalize_Workspace::workspace->v[i]" );
            }

            sfree( workspace->y, "Finalize_Workspace::workspace->y" );
            sfree( workspace->z, "Finalize_Workspace::workspace->z" );
            sfree( workspace->g, "Finalize_Workspace::workspace->g" );
            sfree( workspace->h, "Finalize_Workspace::workspace->h" );
            sfree( workspace->hs, "Finalize_Workspace::workspace->hs" );
            sfree( workspace->hc, "Finalize_Workspace::workspace->hc" );
            sfree( workspace->rn, "Finalize_Workspace::workspace->rn" );
            sfree( workspace->v, "Finalize_Workspace::workspace->v" );

            sfree( workspace->r, "Finalize_Workspace::workspace->r" );
            sfree( workspace->d, "Finalize_Workspace::workspace->d" );
            sfree( workspace->q, "Finalize_Workspace::workspace->q" );
            sfree( workspace->p, "Finalize_Workspace::workspace->p" );
            break;

        /* CG storage */
        case CG_S:
            sfree( workspace->r, "Finalize_Workspace::workspace->r" );
            sfree( workspace->d, "Finalize_Workspace::workspace->d" );
            sfree( workspace->q, "Finalize_Workspace::workspace->q" );
            sfree( workspace->p, "Finalize_Workspace::workspace->p" );
            break;

        case SDM_S:
            sfree( workspace->r, "Finalize_Workspace::workspace->r" );
            sfree( workspace->d, "Finalize_Workspace::workspace->d" );
            sfree( workspace->q, "Finalize_Workspace::workspace->q" );
            break;

        default:
            fprintf( stderr, "Unknown charge method linear solver type. Terminating...\n" );
            exit( INVALID_INPUT );
            break;
    }

    /* SpMV related */
#ifdef _OPENMP
    sfree( workspace->b_local, "Finalize_Workspace::b_local" );
#endif

    /* level scheduling related */
    if ( control->cm_solver_pre_app_type == TRI_SOLVE_LEVEL_SCHED_PA ||
            control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
    {
        sfree( workspace->row_levels_L, "Finalize_Workspace::row_levels_L" );
        sfree( workspace->level_rows_L, "Finalize_Workspace::level_rows_L" );
        sfree( workspace->level_rows_cnt_L, "Finalize_Workspace::level_rows_cnt_L" );
        sfree( workspace->row_levels_U, "Finalize_Workspace::row_levels_U" );
        sfree( workspace->level_rows_U, "Finalize_Workspace::level_rows_U" );
        sfree( workspace->level_rows_cnt_U, "Finalize_Workspace::level_rows_cnt_U" );
        sfree( workspace->top, "Finalize_Workspace::top" );
    }

    /* graph coloring related */
    if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
    {
        sfree( workspace->color, "Finalize_Workspace::workspace->color" );
        sfree( workspace->to_color, "Finalize_Workspace::workspace->to_color" );
        sfree( workspace->conflict, "Finalize_Workspace::workspace->conflict" );
        sfree( workspace->conflict_cnt, "Finalize_Workspace::workspace->conflict_cnt" );
        sfree( workspace->recolor, "Finalize_Workspace::workspace->recolor" );
        sfree( workspace->color_top, "Finalize_Workspace::workspace->color_top" );
        sfree( workspace->permuted_row_col, "Finalize_Workspace::workspace->permuted_row_col" );
        sfree( workspace->permuted_row_col_inv, "Finalize_Workspace::workspace->permuted_row_col_inv" );
        sfree( workspace->y_p, "Finalize_Workspace::workspace->y_p" );
        sfree( workspace->x_p, "Finalize_Workspace::workspace->x_p" );
    }

    /* Jacobi iteration related */
    if ( control->cm_solver_pre_app_type == JACOBI_ITER_PA )
    {
        sfree( workspace->Dinv_L, "Finalize_Workspace::Dinv_L" );
        sfree( workspace->Dinv_U, "Finalize_Workspace::Dinv_U" );
        sfree( workspace->Dinv_b, "Finalize_Workspace::Dinv_b" );
        sfree( workspace->rp, "Finalize_Workspace::rp" );
        sfree( workspace->rp2, "Finalize_Workspace::rp2" );
    }

    /* integrator storage */
    sfree( workspace->a, "Finalize_Workspace::workspace->a" );
    sfree( workspace->f_old, "Finalize_Workspace::workspace->f_old" );
    sfree( workspace->v_const, "Finalize_Workspace::workspace->v_const" );

#ifdef _OPENMP
    sfree( workspace->f_local, "Finalize_Workspace::workspace->f_local" );
#endif

    /* storage for analysis */
    if ( control->molec_anal || control->diffusion_coef )
    {
        sfree( workspace->mark, "Finalize_Workspace::workspace->mark" );
        sfree( workspace->old_mark, "Finalize_Workspace::workspace->old_mark" );
    }

    if ( control->diffusion_coef )
    {
        sfree( workspace->x_old, "Finalize_Workspace::workspace->x_old" );
    }

    sfree( workspace->orig_id, "Finalize_Workspace::workspace->orig_id" );

    /* space for keeping restriction info, if any */
    if ( control->restrict_bonds )
    {
        for ( i = 0; i < system->N; ++i )
        {
            sfree( workspace->restricted_list[i],
                    "Finalize_Workspace::workspace->restricted_list[i]" );
        }

        sfree( workspace->restricted, "Finalize_Workspace::workspace->restricted" );
        sfree( workspace->restricted_list, "Finalize_Workspace::workspace->restricted_list" );
    }

#ifdef TEST_FORCES
    sfree( workspace->dDelta, "Finalize_Workspace::workspace->dDelta" );
    sfree( workspace->f_ele, "Finalize_Workspace::workspace->f_ele" );
    sfree( workspace->f_vdw, "Finalize_Workspace::workspace->f_vdw" );
    sfree( workspace->f_bo, "Finalize_Workspace::workspace->f_bo" );
    sfree( workspace->f_be, "Finalize_Workspace::workspace->f_be" );
    sfree( workspace->f_lp, "Finalize_Workspace::workspace->f_lp" );
    sfree( workspace->f_ov, "Finalize_Workspace::workspace->f_ov" );
    sfree( workspace->f_un, "Finalize_Workspace::workspace->f_un" );
    sfree( workspace->f_ang, "Finalize_Workspace::workspace->f_ang" );
    sfree( workspace->f_coa, "Finalize_Workspace::workspace->f_coa" );
    sfree( workspace->f_pen, "Finalize_Workspace::workspace->f_pen" );
    sfree( workspace->f_hb, "Finalize_Workspace::workspace->f_hb" );
    sfree( workspace->f_tor, "Finalize_Workspace::workspace->f_tor" );
    sfree( workspace->f_con, "Finalize_Workspace::workspace->f_con" );
#endif
}


void Finalize_Lists( control_params *control, reax_list **lists )
{
    Delete_List( TYP_FAR_NEIGHBOR, &(*lists)[FAR_NBRS] );
    if ( control->hb_cut > 0.0 )
    {
        Delete_List( TYP_HBOND, &(*lists)[HBONDS] );
    }
    Delete_List( TYP_BOND, &(*lists)[BONDS] );
    Delete_List( TYP_THREE_BODY, &(*lists)[THREE_BODIES] );

#ifdef TEST_FORCES
    Delete_List( TYP_DDELTA, &(*lists)[DDELTA] );
    Delete_List( TYP_DBO, &(*lists)[DBO] );
#endif
}


void Finalize_Out_Controls( reax_system *system, control_params *control,
        static_storage *workspace, output_controls *out_control )
{
    /* close trajectory file */
    if ( out_control->write_steps > 0 )
    {
        if ( out_control->trj )
        {
            fclose( out_control->trj );
        }
    }

    if ( out_control->energy_update_freq > 0 )
    {
        /* close out file */
        if ( out_control->out )
        {
            fclose( out_control->out );
        }

        /* close potentials file */
        if ( out_control->pot )
        {
            fclose( out_control->pot );
        }

        /* close log file */
        if ( out_control->log )
        {
            fclose( out_control->log );
        }
    }

    /* close pressure file */
    if ( control->ensemble == NPT ||
            control->ensemble == iNPT ||
            control->ensemble == sNPT )
    {
        if ( out_control->prs )
        {
            fclose( out_control->prs );
        }
    }

    /* close molecular analysis file */
    if ( control->molec_anal )
    {
        fclose( out_control->mol );
    }

    /* close electric dipole moment analysis file */
    if ( control->dipole_anal )
    {
        fclose( out_control->dpl );
    }

    /* close diffusion coef analysis file */
    if ( control->diffusion_coef )
    {
        fclose( out_control->drft );
    }


#ifdef TEST_ENERGY
    /* close bond energy file */
    if ( out_control->ebond )
    {
        fclose( out_control->ebond );
    }

    /* close lone-pair energy file */
    if ( out_control->help )
    {
        fclose( out_control->elp );
    }

    /* close overcoordination energy file */
    if ( out_control->eov )
    {
        fclose( out_control->eov );
    }

    /* close undercoordination energy file */
    if ( out_control->eun )
    {
        fclose( out_control->eun );
    }

    /* close angle energy file */
    if ( out_control->eval )
    {
        fclose( out_control->eval );
    }

    /* close penalty energy file */
    if ( out_control->epen )
    {
        fclose( out_control->epen );
    }

    /* close coalition energy file */
    if ( out_control->ecoa )
    {
        fclose( out_control->ecoa );
    }

    /* close hydrogen bond energy file */
    if ( out_control->ehb )
    {
        fclose( out_control->ehb );
    }

    /* close torsion energy file */
    if ( out_control->etor )
    {
        fclose( out_control->etor );
    }

    /* close conjugation energy file */
    if ( out_control->econ )
    {
        fclose( out_control->econ );
    }

    /* close vdWaals energy file */
    if ( out_control->evdw )
    {
        fclose( out_control->evdw );
    }

    /* close coulomb energy file */
    if ( out_control->ecou )
    {
        fclose( out_control->ecou );
    }
#endif


#ifdef TEST_FORCES
    /* close bond orders file */
    if ( out_control->fbo )
    {
        fclose( out_control->fbo );
    }

    /* close bond orders derivatives file */
    if ( out_control->fdbo )
    {
        fclose( out_control->fdbo );
    }

    /* close bond forces file */
    if ( out_control->fbond )
    {
        fclose( out_control->fbond );
    }

    /* close lone-pair forces file */
    if ( out_control->flp )
    {
        fclose( out_control->flp );
    }

    /* close overcoordination forces file */
    if ( out_control->fatom )
    {
        fclose( out_control->fatom );
    }

    /* close angle forces file */
    if ( out_control->f3body )
    {
        fclose( out_control->f3body );
    }

    /* close hydrogen bond forces file */
    if ( out_control->fhb )
    {
        fclose( out_control->fhb );
    }

    /* close torsion forces file */
    if ( out_control->f4body )
    {
        fclose( out_control->f4body );
    }

    /* close nonbonded forces file */
    if ( out_control->fnonb )
    {
        fclose( out_control->fnonb );
    }

    /* close total force file */
    if ( out_control->ftot )
    {
        fclose( out_control->ftot );
    }

    /* close coulomb forces file */
    if ( out_control->ftot2 )
    {
        fclose( out_control->ftot2 );
    }
#endif
}


void Finalize( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace, reax_list **lists,
        output_controls *out_control, const int output_enabled )
{
    if ( control->tabulate )
    {
        Finalize_LR_Lookup_Table( system, control, workspace );
    }

    if ( output_enabled == TRUE )
    {
        Finalize_Out_Controls( system, control, workspace, out_control );
    }

    Finalize_Lists( control, lists );

    Finalize_Workspace( system, control, workspace );

    Finalize_Simulation_Data( system, control, data, out_control );

    Finalize_System( system, control, data );
}