init_md.c

/*----------------------------------------------------------------------
  SerialReax - Reax Force Field Simulator

  Copyright (2010) Purdue University
  Hasan Metin Aktulga, haktulga@cs.purdue.edu
  Joseph Fogarty, jcfogart@mail.usf.edu
  Sagar Pandit, pandit@usf.edu
  Ananth Y Grama, ayg@cs.purdue.edu

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of
  the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  See the GNU General Public License for more details:
  <http://www.gnu.org/licenses/>.
  ----------------------------------------------------------------------*/

#include "init_md.h"

#include "allocate.h"
#include "box.h"
#include "forces.h"
#include "grid.h"
#include "integrate.h"
#include "neighbors.h"
#include "list.h"
#include "lookup.h"
#include "reset_utils.h"
#include "system_props.h"
#include "tool_box.h"
#include "vector.h"


void Generate_Initial_Velocities( reax_system *system, real T )
{
    int i;
    real scale, norm;

    if ( T <= 0.1 )
    {
        for ( i = 0; i < system->N; i++ )
        {
            rvec_MakeZero( system->atoms[i].v );
        }

#if defined(DEBUG)
        fprintf( stderr, "no random velocities...\n" );
#endif
    }
    else
    {
        for ( i = 0; i < system->N; i++ )
        {
            rvec_Random( system->atoms[i].v );

            norm = rvec_Norm_Sqr( system->atoms[i].v );
            scale = SQRT( system->reaxprm.sbp[ system->atoms[i].type ].mass *
                          norm / (3.0 * K_B * T) );

            rvec_Scale( system->atoms[i].v, 1.0 / scale, system->atoms[i].v );

            /*fprintf( stderr, "v = %f %f %f\n",
            system->atoms[i].v[0],system->atoms[i].v[1],system->atoms[i].v[2]);
            fprintf( stderr, "scale = %f\n", scale );
            fprintf( stderr, "v = %f %f %f\n",
            system->atoms[i].v[0],system->atoms[i].v[1],system->atoms[i].v[2]);*/
        }
    }
}


void Init_System( reax_system *system, control_params *control,
        simulation_data *data )
{
    int i;
    rvec dx;

    if ( !control->restart )
    {
        Reset_Atoms( system );
    }

    Compute_Total_Mass( system, data );
    Compute_Center_of_Mass( system, data, stderr );

    /* reposition atoms */
    // just fit the atoms to the periodic box
    if ( control->reposition_atoms == 0 )
    {
        rvec_MakeZero( dx );
    }
    // put the center of mass to the center of the box
    else if ( control->reposition_atoms == 1 )
    {
        rvec_Scale( dx, 0.5, system->box.box_norms );
        rvec_ScaledAdd( dx, -1., data->xcm );
    }
    // put the center of mass to the origin
    else if ( control->reposition_atoms == 2 )
    {
        rvec_Scale( dx, -1., data->xcm );
    }
    else
    {
        fprintf( stderr, "UNKNOWN OPTION: reposition_atoms. Terminating...\n" );
        exit( UNKNOWN_OPTION );
    }

    for ( i = 0; i < system->N; ++i )
    {
        Inc_on_T3( system->atoms[i].x, dx, &(system->box) );
        /*fprintf( stderr, "%6d%2d%8.3f%8.3f%8.3f\n",
          i, system->atoms[i].type,
          system->atoms[i].x[0], system->atoms[i].x[1], system->atoms[i].x[2] );*/
    }

    /* Initialize velocities so that desired init T can be attained */
    if ( !control->restart || (control->restart && control->random_vel) )
    {
        Generate_Initial_Velocities( system, control->T_init );
    }

    Setup_Grid( system );
}


void Init_Simulation_Data( reax_system *system, control_params *control,
        simulation_data *data, output_controls *out_control,
        evolve_function *Evolve )
{

    Reset_Simulation_Data( data );

    if ( !control->restart )
    {
        data->step = 0;
        data->prev_steps = 0;
    }

    switch ( control->ensemble )
    {
    case NVE:
        data->N_f = 3 * system->N;
        *Evolve = Velocity_Verlet_NVE;
        break;


    case NVT:
        data->N_f = 3 * system->N + 1;
        //control->Tau_T = 100 * data->N_f * K_B * control->T_final;
        if ( !control->restart || (control->restart && control->random_vel) )
        {
            data->therm.G_xi = control->Tau_T * (2.0 * data->E_Kin -
                                                 data->N_f * K_B * control->T );
            data->therm.v_xi = data->therm.G_xi * control->dt;
            data->therm.v_xi_old = 0;
            data->therm.xi = 0;
#if defined(DEBUG_FOCUS)
            fprintf( stderr, "init_md: G_xi=%f Tau_T=%f E_kin=%f N_f=%f v_xi=%f\n",
                     data->therm.G_xi, control->Tau_T, data->E_Kin,
                     data->N_f, data->therm.v_xi );
#endif
        }

        *Evolve = Velocity_Verlet_Nose_Hoover_NVT_Klein;
        break;


    case NPT: // Anisotropic NPT
        fprintf( stderr, "THIS OPTION IS NOT YET IMPLEMENTED! TERMINATING...\n" );
        exit( UNKNOWN_OPTION );
        data->N_f = 3 * system->N + 9;
        if ( !control->restart )
        {
            data->therm.G_xi = control->Tau_T * (2.0 * data->E_Kin -
                    data->N_f * K_B * control->T);
            data->therm.v_xi = data->therm.G_xi * control->dt;
            data->iso_bar.eps = 1.0 / 3.0 * LOG( system->box.volume );
            //data->inv_W = 1. / (data->N_f*K_B*control->T*SQR(control->Tau_P));
            //Compute_Pressure( system, data, workspace );
        }
        *Evolve = Velocity_Verlet_Berendsen_Isotropic_NPT;
        break;


    case sNPT: // Semi-Isotropic NPT
        data->N_f = 3 * system->N + 4;
        *Evolve = Velocity_Verlet_Berendsen_SemiIsotropic_NPT;
        break;


    case iNPT: // Isotropic NPT
        data->N_f = 3 * system->N + 2;
        *Evolve = Velocity_Verlet_Berendsen_Isotropic_NPT;
        break;

    case bNVT:
        data->N_f = 3 * system->N + 1;
        *Evolve = Velocity_Verlet_Berendsen_NVT;
        fprintf (stderr, " Initializing Velocity_Verlet_Berendsen_NVT .... \n");
        break;

    default:
        break;
    }

    Compute_Kinetic_Energy( system, data );

    /* init timing info */
    data->timing.start = Get_Time( );
    data->timing.total = data->timing.start;
    data->timing.nbrs = 0;
    data->timing.init_forces = 0;
    data->timing.bonded = 0;
    data->timing.nonb = 0;
    data->timing.cm = ZERO;
    data->timing.cm_sort_mat_rows = ZERO;
    data->timing.cm_solver_pre_comp = ZERO;
    data->timing.cm_solver_pre_app = ZERO;
    data->timing.cm_solver_iters = 0;
    data->timing.cm_solver_spmv = ZERO;
    data->timing.cm_solver_vector_ops = ZERO;
    data->timing.cm_solver_orthog = ZERO;
    data->timing.cm_solver_tri_solve = ZERO;
}


/* Initialize Taper params */
void Init_Taper( control_params *control )
{
    real d1, d7;
    real swa, swa2, swa3;
    real swb, swb2, swb3;

    swa = control->r_low;
    swb = control->r_cut;

    if ( FABS( swa ) > 0.01 )
    {
        fprintf( stderr, "Warning: non-zero value for lower Taper-radius cutoff\n" );
    }

    if ( swb < 0.0 )
    {
        fprintf( stderr, "Negative value for upper Taper-radius cutoff\n" );
        exit( INVALID_INPUT );
    }
    else if ( swb < 5.0 )
    {
        fprintf( stderr, "Warning: low value for upper Taper-radius cutoff:%f\n", swb );
    }

    d1 = swb - swa;
    d7 = POW( d1, 7.0 );
    swa2 = SQR( swa );
    swa3 = CUBE( swa );
    swb2 = SQR( swb );
    swb3 = CUBE( swb );

    control->Tap7 =  20.0 / d7;
    control->Tap6 = -70.0 * (swa + swb) / d7;
    control->Tap5 =  84.0 * (swa2 + 3.0 * swa * swb + swb2) / d7;
    control->Tap4 = -35.0 * (swa3 + 9.0 * swa2 * swb + 9.0 * swa * swb2 + swb3 ) / d7;
    control->Tap3 = 140.0 * (swa3 * swb + 3.0 * swa2 * swb2 + swa * swb3 ) / d7;
    control->Tap2 = -210.0 * (swa3 * swb2 + swa2 * swb3) / d7;
    control->Tap1 = 140.0 * swa3 * swb3 / d7;
    control->Tap0 = (-35.0 * swa3 * swb2 * swb2 + 21.0 * swa2 * swb3 * swb2 +
            7.0 * swa * swb3 * swb3 + swb3 * swb3 * swb ) / d7;
}


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

    /* Allocate space for hydrogen bond list */
    workspace->hbond_index = (int *) malloc( system->N * sizeof( int ) );

    /* bond order related storage  */
    workspace->total_bond_order = (real *) malloc( system->N * sizeof( real ) );
    workspace->Deltap           = (real *) malloc( system->N * sizeof( real ) );
    workspace->Deltap_boc       = (real *) malloc( system->N * sizeof( real ) );
    workspace->dDeltap_self     = (rvec *) malloc( system->N * sizeof( rvec ) );

    workspace->Delta        = (real *) malloc( system->N * sizeof( real ) );
    workspace->Delta_lp         = (real *) malloc( system->N * sizeof( real ) );
    workspace->Delta_lp_temp    = (real *) malloc( system->N * sizeof( real ) );
    workspace->dDelta_lp        = (real *) malloc( system->N * sizeof( real ) );
    workspace->dDelta_lp_temp   = (real *) malloc( system->N * sizeof( real ) );
    workspace->Delta_e          = (real *) malloc( system->N * sizeof( real ) );
    workspace->Delta_boc        = (real *) malloc( system->N * sizeof( real ) );
    workspace->nlp          = (real *) malloc( system->N * sizeof( real ) );
    workspace->nlp_temp         = (real *) malloc( system->N * sizeof( real ) );
    workspace->Clp          = (real *) malloc( system->N * sizeof( real ) );
    workspace->CdDelta          = (real *) malloc( system->N * sizeof( real ) );
    workspace->vlpex        = (real *) malloc( system->N * sizeof( real ) );

    /* charge method storage */
    switch ( control->charge_method )
    {
        case QEQ_CM:
            system->N_cm = system->N;
            break;
        case EE_CM:
            system->N_cm = system->N + 1;
            break;
        case ACKS2_CM:
            system->N_cm = 2 * system->N + 2;
            break;
        default:
            fprintf( stderr, "Unknown charge method type. Terminating...\n" );
            exit( INVALID_INPUT );
            break;
    }

    workspace->H = NULL;
    workspace->H_sp = NULL;
    workspace->L = NULL;
    workspace->H_spar_patt = NULL;
    workspace->H_app_inv = NULL;
    workspace->U = NULL;
    workspace->Hdia_inv = NULL;
    if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
            control->cm_solver_pre_comp_type == ILUT_PAR_PC )
    {
        workspace->droptol  = (real *) calloc( system->N_cm, sizeof( real ) );
    }
    //TODO: check if unused
    //workspace->w        = (real *) calloc( cm_lin_sys_size, sizeof( real ) );
    //TODO: check if unused
    workspace->b        = (real *) calloc( system->N_cm * 2, sizeof( real ) );
    workspace->b_s      = (real *) calloc( system->N_cm, sizeof( real ) );
    workspace->b_t      = (real *) calloc( system->N_cm, sizeof( real ) );
    workspace->b_prc    = (real *) calloc( system->N_cm * 2, sizeof( real ) );
    workspace->b_prm    = (real *) calloc( system->N_cm * 2, sizeof( real ) );
    workspace->s        = (real**) calloc( 5, sizeof( real* ) );
    workspace->t        = (real**) calloc( 5, sizeof( real* ) );
    for ( i = 0; i < 5; ++i )
    {
        workspace->s[i] = (real *) calloc( system->N_cm, sizeof( real ) );
        workspace->t[i] = (real *) calloc( system->N_cm, sizeof( real ) );
    }

    switch ( control->charge_method )
    {
        case QEQ_CM:
            for ( i = 0; i < system->N; ++i )
            {
                workspace->b_s[i] = -system->reaxprm.sbp[ system->atoms[i].type ].chi;
                workspace->b_t[i] = -1.0;

                //TODO: check if unused (redundant)
                workspace->b[i] = -system->reaxprm.sbp[ system->atoms[i].type ].chi;
                workspace->b[i + system->N] = -1.0;
            }
            break;

        case EE_CM:
            for ( i = 0; i < system->N; ++i )
            {
                workspace->b_s[i] = -system->reaxprm.sbp[ system->atoms[i].type ].chi;

                //TODO: check if unused (redundant)
                workspace->b[i] = -system->reaxprm.sbp[ system->atoms[i].type ].chi;
            }

            workspace->b_s[system->N] = control->cm_q_net;
            workspace->b[system->N] = control->cm_q_net;
            break;

        case ACKS2_CM:
            for ( i = 0; i < system->N; ++i )
            {
                workspace->b_s[i] = -system->reaxprm.sbp[ system->atoms[i].type ].chi;

                //TODO: check if unused (redundant)
                workspace->b[i] = -system->reaxprm.sbp[ system->atoms[i].type ].chi;
            }

            workspace->b_s[system->N] = control->cm_q_net;
            workspace->b[system->N] = control->cm_q_net;

            for ( i = system->N + 1; i < system->N_cm; ++i )
            {
                workspace->b_s[i] = 0.0;

                //TODO: check if unused (redundant)
                workspace->b[i] = 0.0;
            }
            break;

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

    switch ( control->cm_solver_type )
    {
        /* GMRES storage */
        case GMRES_S:
        case GMRES_H_S:
            workspace->y  = (real *)  calloc( control->cm_solver_restart + 1, sizeof( real ) );
            workspace->z  = (real *)  calloc( control->cm_solver_restart + 1, sizeof( real ) );
            workspace->g  = (real *)  calloc( control->cm_solver_restart + 1, sizeof( real ) );
            workspace->h  = (real **) calloc( control->cm_solver_restart + 1, sizeof( real*) );
            workspace->hs = (real *)  calloc( control->cm_solver_restart + 1, sizeof( real ) );
            workspace->hc = (real *)  calloc( control->cm_solver_restart + 1, sizeof( real ) );
            workspace->rn = (real **) calloc( control->cm_solver_restart + 1, sizeof( real*) );
            workspace->v  = (real **) calloc( control->cm_solver_restart + 1, sizeof( real*) );

            for ( i = 0; i < control->cm_solver_restart + 1; ++i )
            {
                workspace->h[i]  = (real *) calloc( control->cm_solver_restart + 1, sizeof( real ) );
                workspace->rn[i] = (real *) calloc( system->N_cm * 2, sizeof( real ) );
                workspace->v[i]  = (real *) calloc( system->N_cm, sizeof( real ) );
            }

            workspace->r = (real *) calloc( system->N_cm, sizeof( real ) );
            workspace->d = (real *) calloc( system->N_cm, sizeof( real ) );
            workspace->q = (real *) calloc( system->N_cm, sizeof( real ) );
            workspace->p = (real *) calloc( system->N_cm, sizeof( real ) );
            break;

        /* CG storage */
        case CG_S:
            workspace->r = (real *) calloc( system->N_cm, sizeof( real ) );
            workspace->d = (real *) calloc( system->N_cm, sizeof( real ) );
            workspace->q = (real *) calloc( system->N_cm, sizeof( real ) );
            workspace->p = (real *) calloc( system->N_cm, sizeof( real ) );
            break;

        case SDM_S:
            workspace->r = (real *) calloc( system->N_cm, sizeof( real ) );
            workspace->d = (real *) calloc( system->N_cm, sizeof( real ) );
            workspace->q = (real *) calloc( system->N_cm, sizeof( real ) );
            break;

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

    /* integrator storage */
    workspace->a = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_old = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->v_const = (rvec *) malloc( system->N * sizeof( rvec ) );

#ifdef _OPENMP
    workspace->f_local = (rvec *) malloc( control->num_threads * system->N * sizeof( rvec ) );
#endif

    /* storage for analysis */
    if ( control->molec_anal || control->diffusion_coef )
    {
        workspace->mark = (int *) calloc( system->N, sizeof(int) );
        workspace->old_mark = (int *) calloc( system->N, sizeof(int) );
    }
    else
    {
        workspace->mark = workspace->old_mark = NULL;
    }

    if ( control->diffusion_coef )
    {
        workspace->x_old = (rvec *) calloc( system->N, sizeof( rvec ) );
    }
    else
    {
        workspace->x_old = NULL;
    }

#ifdef TEST_FORCES
    workspace->dDelta = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_ele = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_vdw = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_bo = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_be = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_lp = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_ov = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_un = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_ang = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_coa = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_pen = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_hb = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_tor = (rvec *) malloc( system->N * sizeof( rvec ) );
    workspace->f_con = (rvec *) malloc( system->N * sizeof( rvec ) );
#endif

    workspace->realloc.num_far = -1;
    workspace->realloc.Htop = -1;
    workspace->realloc.hbonds = -1;
    workspace->realloc.bonds = -1;
    workspace->realloc.num_3body = -1;
    workspace->realloc.gcell_atoms = -1;

    Reset_Workspace( system, workspace );

    /* Initialize Taper function */
    Init_Taper( control );
}


void Init_Lists( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace,
        reax_list **lists, output_controls *out_control )
{
    int i, num_nbrs, num_bonds, num_3body, Htop, max_nnz;
    int *hb_top, *bond_top;
#if defined(DEBUG_FOCUS)
    int num_hbonds;
#endif

    num_nbrs = Estimate_NumNeighbors( system, control, workspace, lists );

    Make_List( system->N, num_nbrs, TYP_FAR_NEIGHBOR, (*lists) + FAR_NBRS );

#if defined(DEBUG_FOCUS)
    fprintf( stderr, "memory allocated: far_nbrs = %ldMB\n",
             num_nbrs * sizeof(far_neighbor_data) / (1024 * 1024) );
#endif

    Generate_Neighbor_Lists(system, control, data, workspace, lists, out_control);
    Htop = 0;
    hb_top = (int*) calloc( system->N, sizeof(int) );
    bond_top = (int*) calloc( system->N, sizeof(int) );
    num_3body = 0;
    Estimate_Storage_Sizes( system, control, lists, &Htop,
            hb_top, bond_top, &num_3body );
    num_3body = MAX( num_3body, MIN_BONDS );

    switch ( control->charge_method )
    {
        case QEQ_CM:
            max_nnz = Htop;
            break;
        case EE_CM:
            max_nnz = Htop + system->N_cm;
            break;
        case ACKS2_CM:
            max_nnz = 2 * Htop + 3 * system->N + 2;
            break;
        default:
            max_nnz = Htop;
            break;
    }

    if ( Allocate_Matrix( &(workspace->H), system->N_cm, max_nnz ) == FAILURE )
    {
        fprintf( stderr, "Not enough space for init matrices. Terminating...\n" );
        exit( INSUFFICIENT_MEMORY );
    }
    /* TODO: better estimate for H_sp?
     *   If so, need to refactor Estimate_Storage_Sizes
     *   to use various cut-off distances as parameters
     *   (non-bonded, hydrogen, 3body, etc.) */
    if ( Allocate_Matrix( &(workspace->H_sp), system->N_cm, max_nnz ) == FAILURE )
    {
        fprintf( stderr, "Not enough space for init matrices. Terminating...\n" );
        exit( INSUFFICIENT_MEMORY );
    }

#if defined(DEBUG_FOCUS)
    fprintf( stderr, "estimated storage - Htop: %d\n", Htop );
    fprintf( stderr, "memory allocated: H = %ldMB\n",
            Htop * sizeof(sparse_matrix_entry) / (1024 * 1024) );
#endif

    workspace->num_H = 0;
    if ( control->hb_cut > 0.0 )
    {
        /* init H indexes */
        for ( i = 0; i < system->N; ++i )
        {
            // H atom
            if ( system->reaxprm.sbp[ system->atoms[i].type ].p_hbond == 1 )
            {
                workspace->hbond_index[i] = workspace->num_H++;
            }
            else
            {
                workspace->hbond_index[i] = -1;
            }
        }

        if ( workspace->num_H == 0 )
        {
            control->hb_cut = 0.0;
        }
        else
        {
            Allocate_HBond_List( system->N, workspace->num_H, workspace->hbond_index,
                    hb_top, (*lists) + HBONDS );
        }

#if defined(DEBUG_FOCUS)
        num_hbonds = hb_top[system->N - 1];
        fprintf( stderr, "estimated storage - num_hbonds: %d\n", num_hbonds );
        fprintf( stderr, "memory allocated: hbonds = %ldMB\n",
                 num_hbonds * sizeof(hbond_data) / (1024 * 1024) );
#endif
    }

    /* bonds list */
    Allocate_Bond_List( system->N, bond_top, (*lists) + BONDS );
    num_bonds = bond_top[system->N - 1];

#if defined(DEBUG_FOCUS)
    fprintf( stderr, "estimated storage - num_bonds: %d\n", num_bonds );
    fprintf( stderr, "memory allocated: bonds = %ldMB\n",
             num_bonds * sizeof(bond_data) / (1024 * 1024) );
#endif

    /* 3bodies list */
    Make_List( num_bonds, num_3body, TYP_THREE_BODY, (*lists) + THREE_BODIES );

#if defined(DEBUG_FOCUS)
    fprintf( stderr, "estimated storage - num_3body: %d\n", num_3body );
    fprintf( stderr, "memory allocated: 3-body = %ldMB\n",
             num_3body * sizeof(three_body_interaction_data) / (1024 * 1024) );
#endif

#ifdef TEST_FORCES
    Make_List( system->N, num_bonds * 8, TYP_DDELTA, (*lists) + DDELTA );

    Make_List( num_bonds, num_bonds * MAX_BONDS * 3, TYP_DBO, (*lists) + DBO );
#endif

    sfree( hb_top, "Init_Lists::hb_top" );
    sfree( bond_top, "Init_Lists::bond_top" );
}


void Init_Out_Controls( reax_system *system, control_params *control,
        static_storage *workspace, output_controls *out_control )
{
#define TEMP_SIZE (1000)
    char temp[TEMP_SIZE];

    /* Init trajectory file */
    if ( out_control->write_steps > 0 )
    {
        strcpy( temp, control->sim_name );
        strcat( temp, ".trj" );
        out_control->trj = fopen( temp, "w" );
        out_control->write_header( system, control, workspace, out_control );
    }

    if ( out_control->energy_update_freq > 0 )
    {
        /* Init out file */
        strcpy( temp, control->sim_name );
        strcat( temp, ".out" );
        out_control->out = fopen( temp, "w" );
        fprintf( out_control->out, "%-6s%16s%16s%16s%11s%11s%13s%13s%13s\n",
                 "step", "total energy", "poten. energy", "kin. energy",
                 "temp.", "target", "volume", "press.", "target" );
        fflush( out_control->out );

        /* Init potentials file */
        strcpy( temp, control->sim_name );
        strcat( temp, ".pot" );
        out_control->pot = fopen( temp, "w" );
        fprintf( out_control->pot,
                 "%-6s%13s%13s%13s%13s%13s%13s%13s%13s%13s%13s%13s\n",
                 "step", "ebond", "eatom", "elp", "eang", "ecoa", "ehb",
                 "etor", "econj", "evdw", "ecoul", "epol" );
        fflush( out_control->pot );

        /* Init log file */
        strcpy( temp, control->sim_name );
        strcat( temp, ".log" );
        out_control->log = fopen( temp, "w" );
        fprintf( out_control->log, "%-6s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n",
                 "step", "total", "neighbors", "init", "bonded",
                 "nonbonded", "CM", "CM Sort", "S iters", "Pre Comp", "Pre App",
                 "S spmv", "S vec ops", "S orthog", "S tsolve" );
    }

    /* Init pressure file */
    if ( control->ensemble == NPT ||
            control->ensemble == iNPT ||
            control->ensemble == sNPT )
    {
        strcpy( temp, control->sim_name );
        strcat( temp, ".prs" );
        out_control->prs = fopen( temp, "w" );
        fprintf( out_control->prs, "%-6s%13s%13s%13s%13s%13s%13s%13s%13s\n",
                 "step", "norm_x", "norm_y", "norm_z",
                 "press_x", "press_y", "press_z", "target_p", "volume" );
        fflush( out_control->prs );
    }

    /* Init molecular analysis file */
    if ( control->molec_anal )
    {
        snprintf( temp, TEMP_SIZE + 4, "%s.mol", control->sim_name );
        out_control->mol = fopen( temp, "w" );
        if ( control->num_ignored )
        {
            snprintf( temp, TEMP_SIZE + 4, "%s.ign", control->sim_name );
            out_control->ign = fopen( temp, "w" );
        }
    }

    /* Init electric dipole moment analysis file */
    if ( control->dipole_anal )
    {
        strcpy( temp, control->sim_name );
        strcat( temp, ".dpl" );
        out_control->dpl = fopen( temp, "w" );
        fprintf( out_control->dpl,
                 "Step      Molecule Count  Avg. Dipole Moment Norm\n" );
        fflush( out_control->dpl );
    }

    /* Init diffusion coef analysis file */
    if ( control->diffusion_coef )
    {
        strcpy( temp, control->sim_name );
        strcat( temp, ".drft" );
        out_control->drft = fopen( temp, "w" );
        fprintf( out_control->drft, "Step     Type Count   Avg Squared Disp\n" );
        fflush( out_control->drft );
    }


#ifdef TEST_ENERGY
    /* open bond energy file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".ebond" );
    out_control->ebond = fopen( temp, "w" );

    /* open lone-pair energy file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".elp" );
    out_control->elp = fopen( temp, "w" );

    /* open overcoordination energy file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".eov" );
    out_control->eov = fopen( temp, "w" );

    /* open undercoordination energy file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".eun" );
    out_control->eun = fopen( temp, "w" );

    /* open angle energy file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".eval" );
    out_control->eval = fopen( temp, "w" );

    /* open penalty energy file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".epen" );
    out_control->epen = fopen( temp, "w" );

    /* open coalition energy file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".ecoa" );
    out_control->ecoa = fopen( temp, "w" );

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

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

    /* open conjugation energy file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".econ" );
    out_control->econ = fopen( temp, "w" );

    /* open vdWaals energy file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".evdw" );
    out_control->evdw = fopen( temp, "w" );

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


#ifdef TEST_FORCES
    /* open bond orders file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".fbo" );
    out_control->fbo = fopen( temp, "w" );

    /* open bond orders derivatives file */
    strcpy( temp, control->sim_name );
    strcat( temp, ".fdbo" );
    out_control->fdbo = fopen( temp, "w" );

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

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

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

    /* open angle forces file */
    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,
        const int output_enabled )
{
#if defined(DEBUG)
    real start, end;
#endif

#ifdef _OPENMP
    #pragma omp parallel default(shared)
    {
        #pragma omp single
        control->num_threads = omp_get_num_threads( );
    }
#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 );

#ifdef TEST_FORCES
    Init_Force_Test_Functions( );
#endif

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

        Make_LR_Lookup_Table( system, control );

#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_spar_patt );
        Deallocate_Matrix( workspace->H_app_inv );
    }

    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;
    }

    /* 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 );
    }

    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 );
}