forces.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 "forces.h"

#include "box.h"
#include "bond_orders.h"
#include "bonds.h"
#include "charges.h"
#include "hydrogen_bonds.h"
#if defined(TEST_FORCES)
  #include "io_tools.h"
#endif
#include "list.h"
#include "multi_body.h"
#include "nonbonded.h"
#include "system_props.h"
#include "torsion_angles.h"
#include "tool_box.h"
#include "valence_angles.h"
#include "vector.h"


typedef enum
{
    DIAGONAL = 0,
    OFF_DIAGONAL = 1,
} MATRIX_ENTRY_POSITION;


#if defined(TEST_FORCES)
static int compare_bonds( const void *p1, const void *p2 )
{
    return ((bond_data *)p1)->nbr - ((bond_data *)p2)->nbr;
}
#endif


void Init_Bonded_Force_Functions( control_params *control )
{
    control->intr_funcs[0] = &BO;
    control->intr_funcs[1] = &Bonds;
    control->intr_funcs[2] = &Atom_Energy;
    control->intr_funcs[3] = &Valence_Angles;
    control->intr_funcs[4] = &Torsion_Angles;
    if ( control->hbond_cut > 0.0 )
    {
        control->intr_funcs[5] = &Hydrogen_Bonds;
    }
    else
    {
        control->intr_funcs[5] = NULL;
    }
    control->intr_funcs[6] = NULL;
    control->intr_funcs[7] = NULL;
    control->intr_funcs[8] = NULL;
    control->intr_funcs[9] = NULL;
}


static void Compute_Bonded_Forces( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace, reax_list **lists,
        output_controls *out_control )
{
    int i;

#if defined(TEST_ENERGY)
    /* Mark beginning of a new timestep in each energy file */
    fprintf( out_control->ebond, "step: %d\n%6s%6s%12s%12s%12s\n",
             data->step, "atom1", "atom2", "bo", "ebond", "total" );
    fprintf( out_control->elp, "step: %d\n%6s%12s%12s%12s\n",
             data->step, "atom", "nlp", "elp", "total" );
    fprintf( out_control->eov, "step: %d\n%6s%12s%12s\n",
             data->step, "atom", "eov", "total" );
    fprintf( out_control->eun, "step: %d\n%6s%12s%12s\n",
             data->step, "atom", "eun", "total" );
    fprintf( out_control->eval, "step: %d\n%6s%6s%6s%12s%12s%12s%12s%12s%12s\n",
             data->step, "atom1", "atom2", "atom3",
             "angle", "bo(12)", "bo(23)", "eval", "epen", "total" );
    fprintf( out_control->epen, "step: %d\n%6s%6s%6s%12s%12s%12s%12s%12s\n",
             data->step, "atom1", "atom2", "atom3",
             "angle", "bo(12)", "bo(23)", "epen", "total" );
    fprintf( out_control->ecoa, "step: %d\n%6s%6s%6s%12s%12s%12s%12s%12s\n",
             data->step, "atom1", "atom2", "atom3",
             "angle", "bo(12)", "bo(23)", "ecoa", "total" );
    fprintf( out_control->ehb,  "step: %d\n%6s%6s%6s%12s%12s%12s%12s%12s\n",
             data->step, "atom1", "atom2", "atom3",
             "r(23)", "angle", "bo(12)", "ehb", "total" );
    fprintf( out_control->etor, "step: %d\n%6s%6s%6s%6s%12s%12s%12s%12s\n",
             data->step, "atom1", "atom2", "atom3", "atom4",
             "phi", "bo(23)", "etor", "total" );
    fprintf( out_control->econ, "step:%d\n%6s%6s%6s%6s%12s%12s%12s%12s%12s%12s\n",
             data->step, "atom1", "atom2", "atom3", "atom4",
             "phi", "bo(12)", "bo(23)", "bo(34)", "econ", "total" );
#endif

    /* function calls for bonded interactions */
    for ( i = 0; i < NUM_INTRS; i++ )
    {
        if ( control->intr_funcs[i] != NULL )
        {
            (control->intr_funcs[i])( system, control, data, workspace,
                    lists, out_control );
        }
    }

#if defined(TEST_FORCES)
    /* function calls for printing bonded interactions */
    for ( i = 0; i < NUM_INTRS; i++ )
    {
        if ( control->print_intr_funcs[i] != NULL )
        {
            (control->print_intr_funcs[i])( system, control, data, workspace,
                    lists, out_control );
        }
    }
#endif
}


static void Compute_NonBonded_Forces( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace,
        reax_list** lists, output_controls *out_control )
{
    real t_start, t_elapsed;

#if defined(TEST_ENERGY)
    fprintf( out_control->evdw, "step: %d\n%6s%6s%12s%12s%12s\n",
             data->step, "atom1", "atom2", "r12", "evdw", "total" );
    fprintf( out_control->ecou, "step: %d\n%6s%6s%12s%12s%12s%12s%12s\n",
             data->step, "atom1", "atom2", "r12", "q1", "q2", "ecou", "total" );
#endif

    t_start = Get_Time( );
    Compute_Charges( system, control, data, workspace, out_control );
    t_elapsed = Get_Timing_Info( t_start );
    data->timing.cm += t_elapsed;
    
    if ( control->cm_solver_pre_comp_refactor == -1 )
    {
        if ( data->step <= 4 || is_refactoring_step( control, data ) )
        {
            if ( is_refactoring_step( control, data ) )
            {
                data->timing.cm_last_pre_comp = data->timing.cm_solver_pre_comp;
            }

            data->timing.cm_optimum = data->timing.cm_solver_pre_app
                + data->timing.cm_solver_spmv
                + data->timing.cm_solver_vector_ops
                + data->timing.cm_solver_orthog
                + data->timing.cm_solver_tri_solve;
            data->timing.cm_total_loss = ZERO;
        }
        else
        {
            data->timing.cm_total_loss += data->timing.cm_solver_pre_app
                + data->timing.cm_solver_spmv
                + data->timing.cm_solver_vector_ops
                + data->timing.cm_solver_orthog
                + data->timing.cm_solver_tri_solve
                - data->timing.cm_optimum;
        }
    }

    if ( control->tabulate <= 0 )
    {
        vdW_Coulomb_Energy( system, control, data, workspace, lists, out_control );
    }
    else
    {
        Tabulated_vdW_Coulomb_Energy( system, control, data, workspace,
                lists, out_control );
    }

#if defined(TEST_FORCES)
    Print_vdW_Coulomb_Forces( system, control, data, workspace,
            lists, out_control );
#endif
}


/* This version of Compute_Total_Force computes forces from coefficients
   accumulated by all interaction functions. Saves enormous time & space! */
static void Compute_Total_Force( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace, reax_list **lists )
{
    int i;
    reax_list *bonds;

    bonds = lists[BONDS];

#ifdef _OPENMP
    #pragma omp parallel default(shared)
#endif
    {
        int pj;
#ifdef _OPENMP
        int j;
#endif

#ifdef _OPENMP
        #pragma omp for schedule(static)
#endif
        for ( i = 0; i < system->N; ++i )
        {
            for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
            {
                if ( i < bonds->bond_list[pj].nbr )
                {
                    if ( control->ensemble == NVE || control->ensemble == nhNVT
                            || control->ensemble == bNVT )
                    {
                        Add_dBond_to_Forces( i, pj, system, data, workspace, lists );
                    }
                    else
                    {
                        Add_dBond_to_Forces_NPT( i, pj, system, data, workspace, lists );
                    }
                }
            }
        }

#ifdef _OPENMP
        /* reduction (sum) on thread-local force vectors */
        #pragma omp for schedule(static)
        for ( i = 0; i < system->N; ++i )
        {
            for ( j = 0; j < control->num_threads; ++j )
            {
                rvec_Add( system->atoms[i].f, workspace->f_local[j * system->N + i] );
            }
        }
#endif
    }
}


static void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int n,
        int Hmax, int Htop, int num_bonds, int num_hbonds )
{
    int i, flag;
    reax_list *bonds, *hbonds;

    bonds = lists[BONDS];
    hbonds = lists[HBONDS];

    /* far neighbors */
    if ( Htop > Hmax * DANGER_ZONE )
    {
        workspace->realloc.Htop = Htop;
        if ( Htop > Hmax )
        {
            fprintf( stderr,
                     "[ERROR] step%d - ran out of space on H matrix: Htop=%d, max = %d",
                     step, Htop, Hmax );
            exit( INSUFFICIENT_MEMORY );
        }
    }

    /* bond list */
    flag = -1;
    workspace->realloc.num_bonds = num_bonds;
    for ( i = 0; i < n - 1; ++i )
    {
        if ( End_Index(i, bonds) >= Start_Index(i + 1, bonds) - 2 )
        {
            workspace->realloc.bonds = 1;
            if ( End_Index(i, bonds) > Start_Index(i + 1, bonds) )
            {
                flag = i;
            }
        }
    }

    if ( flag > -1 )
    {
        fprintf( stderr, "[ERROR] step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
                 step, flag, End_Index(flag, bonds), Start_Index(flag + 1, bonds) );
        exit( INSUFFICIENT_MEMORY );
    }

    if ( End_Index(i, bonds) >= bonds->total_intrs - 2 )
    {
        workspace->realloc.bonds = 1;

        if ( End_Index(i, bonds) > bonds->total_intrs )
        {
            fprintf( stderr, "[ERROR] step%d-bondchk failed: i=%d end(i)=%d bond_end=%d\n",
                     step, flag, End_Index(i, bonds), bonds->total_intrs );
            exit( INSUFFICIENT_MEMORY );
        }
    }


    /* hbonds list */
    if ( workspace->num_H > 0 )
    {
        flag = -1;
        workspace->realloc.num_hbonds = num_hbonds;
        for ( i = 0; i < workspace->num_H - 1; ++i )
        {
            if ( Num_Entries(i, hbonds) >=
                    (Start_Index(i + 1, hbonds) - Start_Index(i, hbonds)) * DANGER_ZONE )
            {
                workspace->realloc.hbonds = 1;
                if ( End_Index(i, hbonds) > Start_Index(i + 1, hbonds) )
                {
                    flag = i;
                }
            }
        }

        if ( flag > -1 )
        {
            fprintf( stderr, "[ERROR] step%d-hbondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
                     step, flag, End_Index(flag, hbonds), Start_Index(flag + 1, hbonds) );
            exit( INSUFFICIENT_MEMORY );
        }

        if ( Num_Entries(i, hbonds) >=
                (hbonds->total_intrs - Start_Index(i, hbonds)) * DANGER_ZONE )
        {
            workspace->realloc.hbonds = 1;

            if ( End_Index(i, hbonds) > hbonds->total_intrs )
            {
                fprintf( stderr, "[ERROR] step%d-hbondchk failed: i=%d end(i)=%d hbondend=%d\n",
                         step, flag, End_Index(i, hbonds), hbonds->total_intrs );
                exit( INSUFFICIENT_MEMORY );
            }
        }
    }
}


static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
        control_params *control, LR_lookup_table **LR, int i, int j,
        real r_ij, MATRIX_ENTRY_POSITION pos )
{
    int r;
    real base, dif, val, ret;
    LR_lookup_table *t;

    ret = 0.0;

    switch ( control->charge_method )
    {
    case QEQ_CM:
    //TODO: tabulate other portions of matrices for EE, ACKS2?
    case EE_CM:
    case ACKS2_CM:
        switch ( pos )
        {
            case OFF_DIAGONAL:
                t = &LR[MIN( system->atoms[i].type, system->atoms[j].type )]
                       [MAX( system->atoms[i].type, system->atoms[j].type )];

                /* cubic spline interpolation */
                r = (int)(r_ij * t->inv_dx);
                if ( r == 0 ) 
                {
                    ++r;
                }
                base = (real)(r + 1) * t->dx;
                dif = r_ij - base;
                val = ((t->ele[r].d * dif + t->ele[r].c) * dif + t->ele[r].b)
                    * dif + t->ele[r].a;
                val *= EV_to_KCALpMOL / C_ELE;

                ret = ((i == j) ? 0.5 : 1.0) * val;
            break;

            case DIAGONAL:
                ret = system->reax_param.sbp[system->atoms[i].type].eta;
            break;

            default:
                fprintf( stderr, "[ERROR] Invalid matrix position. Terminating...\n" );
                exit( INVALID_INPUT );
            break;
        }
        break;

    default:
        fprintf( stderr, "[ERROR] Invalid charge method. Terminating...\n" );
        exit( INVALID_INPUT );
        break;
    }

    return ret;
}


static inline real Init_Charge_Matrix_Entry( reax_system *system,
        control_params *control, static_storage *workspace,
        int i, int j, real r_ij, MATRIX_ENTRY_POSITION pos )
{
    real Tap, dr3gamij_1, dr3gamij_3, ret;

    ret = 0.0;

    switch ( control->charge_method )
    {
    case QEQ_CM:
    case EE_CM:
        switch ( pos )
        {
            case OFF_DIAGONAL:
                Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
                Tap = Tap * r_ij + workspace->Tap[5];
                Tap = Tap * r_ij + workspace->Tap[4];
                Tap = Tap * r_ij + workspace->Tap[3];
                Tap = Tap * r_ij + workspace->Tap[2];
                Tap = Tap * r_ij + workspace->Tap[1];
                Tap = Tap * r_ij + workspace->Tap[0];

                /* shielding */
                dr3gamij_1 = r_ij * r_ij * r_ij
                        + POW( system->reax_param.tbp[system->atoms[i].type][system->atoms[j].type].gamma, -3.0 );
                dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );

                /* i == j: periodic self-interaction term
                 * i != j: general interaction term */
                ret = ((i == j) ? 0.5 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
            break;

            case DIAGONAL:
                ret = system->reax_param.sbp[system->atoms[i].type].eta;
            break;

            default:
                fprintf( stderr, "[ERROR] Invalid matrix position. Terminating...\n" );
                exit( INVALID_INPUT );
            break;
        }
        break;

    case ACKS2_CM:
        switch ( pos )
        {
            case OFF_DIAGONAL:
                Tap = workspace->Tap[7] * r_ij
                    + workspace->Tap[6];
                Tap = Tap * r_ij + workspace->Tap[5];
                Tap = Tap * r_ij + workspace->Tap[4];
                Tap = Tap * r_ij + workspace->Tap[3];
                Tap = Tap * r_ij + workspace->Tap[2];
                Tap = Tap * r_ij + workspace->Tap[1];
                Tap = Tap * r_ij + workspace->Tap[0];

                /* shielding */
                dr3gamij_1 = r_ij * r_ij * r_ij
                        + POW( system->reax_param.tbp[system->atoms[i].type][system->atoms[j].type].gamma, -3.0 );
                dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );

                /* i == j: periodic self-interaction term
                 * i != j: general interaction term */
                ret = ((i == j) ? 0.0 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
            break;

            case DIAGONAL:
                /* parameters in electron-volts */
                ret = system->reax_param.sbp[system->atoms[i].type].eta;
            break;

            default:
                fprintf( stderr, "[ERROR] Invalid matrix position. Terminating...\n" );
                exit( INVALID_INPUT );
            break;
        }
        break;

    default:
        fprintf( stderr, "[ERROR] Invalid charge method. Terminating...\n" );
        exit( INVALID_INPUT );
        break;
    }

    return ret;
}


static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
        control_params *control, reax_list *far_nbrs,
        sparse_matrix * H, sparse_matrix * H_sp,
        int * Htop, int * H_sp_top )
{
    int i, j, pj, target, val_flag;
    real d, xcut, bond_softness, * X_diag;

    switch ( control->charge_method )
    {
        case QEQ_CM:
            break;

        case EE_CM:
            H->start[system->N_cm - 1] = *Htop;
            H_sp->start[system->N_cm - 1] = *H_sp_top;

            for ( i = 0; i < system->N_cm - 1; ++i )
            {
                H->j[*Htop] = i;
                H->val[*Htop] = 1.0;
                *Htop = *Htop + 1;

                H_sp->j[*H_sp_top] = i;
                H_sp->val[*H_sp_top] = 1.0;
                *H_sp_top = *H_sp_top + 1;
            }

            H->j[*Htop] = system->N_cm - 1;
            H->val[*Htop] = 0.0;
            *Htop = *Htop + 1;

            H_sp->j[*H_sp_top] = system->N_cm - 1;
            H_sp->val[*H_sp_top] = 0.0;
            *H_sp_top = *H_sp_top + 1;
            break;

        case ACKS2_CM:
            X_diag = smalloc( sizeof(real) * system->N,
                    "Init_Charge_Matrix_Remaining_Entries::X_diag" );

            for ( i = 0; i < system->N; ++i )
            {
                X_diag[i] = 0.0;
            }

            for ( i = 0; i < system->N; ++i )
            {
                H->start[system->N + i] = *Htop;
                H_sp->start[system->N + i] = *H_sp_top;

                /* constraint on ref. value for kinetic energy potential */
                H->j[*Htop] = i;
                H->val[*Htop] = 1.0;
                *Htop = *Htop + 1;

                H_sp->j[*H_sp_top] = i;
                H_sp->val[*H_sp_top] = 1.0;
                *H_sp_top = *H_sp_top + 1;

                /* kinetic energy terms */
                for ( pj = Start_Index(i, far_nbrs); pj < End_Index(i, far_nbrs); ++pj )
                {
                    /* exclude self-periodic images of atoms for
                     * kinetic energy term because the effective
                     * potential is the same on an atom and its periodic image */
                    if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
                    {
                        j = far_nbrs->far_nbr_list[pj].nbr;

                        xcut = 0.5 * ( system->reax_param.sbp[ system->atoms[i].type ].b_s_acks2
                                + system->reax_param.sbp[ system->atoms[j].type ].b_s_acks2 );

                        if ( far_nbrs->far_nbr_list[pj].d < xcut )
                        {
                            d = far_nbrs->far_nbr_list[pj].d / xcut;
                            bond_softness = system->reax_param.gp.l[34] * POW( d, 3.0 )
                                * POW( 1.0 - d, 6.0 );

                            if ( bond_softness > 0.0 )
                            {
                                val_flag = FALSE;

                                for ( target = H->start[system->N + i]; target < *Htop; ++target )
                                {
                                    if ( H->j[target] == system->N + j )
                                    {
                                        H->val[target] += bond_softness;
                                        val_flag = TRUE;
                                        break;
                                    }
                                }

                                if ( val_flag == FALSE )
                                {
                                    H->j[*Htop] = system->N + j;
                                    H->val[*Htop] = bond_softness;
                                    ++(*Htop);
                                }

                                val_flag = FALSE;

                                for ( target = H_sp->start[system->N + i]; target < *H_sp_top; ++target )
                                {
                                    if ( H_sp->j[target] == system->N + j )
                                    {
                                        H_sp->val[target] += bond_softness;
                                        val_flag = TRUE;
                                        break;
                                    }
                                }

                                if ( val_flag == FALSE )
                                {
                                    H_sp->j[*H_sp_top] = system->N + j;
                                    H_sp->val[*H_sp_top] = bond_softness;
                                    ++(*H_sp_top);
                                }

                                X_diag[i] -= bond_softness;
                                X_diag[j] -= bond_softness;
                            }
                        }
                    }
                }

                /* placeholders for diagonal entries, to be replaced below */
                H->j[*Htop] = system->N + i;
                H->val[*Htop] = 0.0;
                *Htop = *Htop + 1;

                H_sp->j[*H_sp_top] = system->N + i;
                H_sp->val[*H_sp_top] = 0.0;
                *H_sp_top = *H_sp_top + 1;
            }

            /* second to last row */
            H->start[system->N_cm - 2] = *Htop;
            H_sp->start[system->N_cm - 2] = *H_sp_top;

            /* place accumulated diagonal entries (needed second to last row marker above before this code) */
            for ( i = system->N; i < 2 * system->N; ++i )
            {
                for ( pj = H->start[i]; pj < H->start[i + 1]; ++pj )
                {
                    if ( H->j[pj] == i )
                    {
                        H->val[pj] = X_diag[i - system->N];
                        break;
                    }
                }

                for ( pj = H_sp->start[i]; pj < H_sp->start[i + 1]; ++pj )
                {
                    if ( H_sp->j[pj] == i )
                    {
                        H_sp->val[pj] = X_diag[i - system->N];
                        break;
                    }
                }
            }

            /* coupling with the kinetic energy potential */
            for ( i = 0; i < system->N; ++i )
            {
                H->j[*Htop] = system->N + i;
                H->val[*Htop] = 1.0;
                *Htop = *Htop + 1;

                H_sp->j[*H_sp_top] = system->N + i;
                H_sp->val[*H_sp_top] = 1.0;
                *H_sp_top = *H_sp_top + 1;
            }

            /* explicitly store zero on diagonal */
            H->j[*Htop] = system->N_cm - 2;
            H->val[*Htop] = 0.0;
            *Htop = *Htop + 1;

            H_sp->j[*H_sp_top] = system->N_cm - 2;
            H_sp->val[*H_sp_top] = 0.0;
            *H_sp_top = *H_sp_top + 1;

            /* last row */
            H->start[system->N_cm - 1] = *Htop;
            H_sp->start[system->N_cm - 1] = *H_sp_top;

            for ( i = 0; i < system->N; ++i )
            {
                H->j[*Htop] = i;
                H->val[*Htop] = 1.0;
                *Htop = *Htop + 1;

                H_sp->j[*H_sp_top] = i;
                H_sp->val[*H_sp_top] = 1.0;
                *H_sp_top = *H_sp_top + 1;
            }

            /* explicitly store zero on diagonal */
            H->j[*Htop] = system->N_cm - 1;
            H->val[*Htop] = 0.0;
            *Htop = *Htop + 1;

            H_sp->j[*H_sp_top] = system->N_cm - 1;
            H_sp->val[*H_sp_top] = 0.0;
            *H_sp_top = *H_sp_top + 1;

            sfree( X_diag, "Init_Charge_Matrix_Remaining_Entries::X_diag" );
            break;

        default:
            break;
    }
}


static void Init_Forces( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace,
        reax_list **lists, output_controls *out_control )
{
    int i, j, pj, target;
    int start_i, end_i;
    int type_i, type_j;
    int Htop, H_sp_top, btop_i, btop_j, num_bonds, num_hbonds;
    int ihb, jhb, ihb_top, jhb_top;
    int flag, flag_sp, val_flag, renbr;
    real r_ij, r2, val;
    real C12, C34, C56;
    real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
    real BO, BO_s, BO_pi, BO_pi2;
    sparse_matrix *H, *H_sp;
    reax_list *far_nbrs, *bonds, *hbonds;
    single_body_parameters *sbp_i, *sbp_j;
    two_body_parameters *twbp;
    far_neighbor_data *nbr_pj;
    reax_atom *atom_i, *atom_j;
    bond_data *ibond, *jbond;
    bond_order_data *bo_ij, *bo_ji;

    far_nbrs = lists[FAR_NBRS];
    bonds = lists[BONDS];
    hbonds = lists[HBONDS];
    H = workspace->H;
    H_sp = workspace->H_sp;
    Htop = 0;
    H_sp_top = 0;
    num_bonds = 0;
    num_hbonds = 0;
    btop_i = 0;
    btop_j = 0;
    renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;

    for ( i = 0; i < system->N; ++i )
    {
        atom_i = &system->atoms[i];
        type_i = atom_i->type;
        start_i = Start_Index( i, far_nbrs );
        end_i = End_Index( i, far_nbrs );
        H->start[i] = Htop;
        H_sp->start[i] = H_sp_top;
        btop_i = End_Index( i, bonds );
        sbp_i = &system->reax_param.sbp[type_i];

        if ( control->hbond_cut > 0.0 )
        {
            ihb = sbp_i->p_hbond;

            if ( ihb == H_ATOM )
            {
                ihb_top = End_Index( workspace->hbond_index[i], hbonds );
            }
            else
            {
                ihb_top = -1;
            }
        }
        else
        {
            ihb = NON_H_BONDING_ATOM;
            ihb_top = -1;
        }

        for ( pj = start_i; pj < end_i; ++pj )
        {
            nbr_pj = &far_nbrs->far_nbr_list[pj];
            j = nbr_pj->nbr;
            flag = FALSE;
            flag_sp = FALSE;

            /* check if reneighboring step --
             * atomic distances just computed via
             * Verlet list, so use current distances */
            if ( renbr == TRUE )
            {
                if ( nbr_pj->d <= control->nonb_cut )
                {
                    flag = TRUE;

                    if ( nbr_pj->d <= control->nonb_sp_cut )
                    {
                        flag_sp = TRUE;
                    }
                }
            }
            /* update atomic distances */
            else
            {
                atom_j = &system->atoms[j];
                nbr_pj->d = Compute_Atom_Distance( &system->box,
                        atom_i->x, atom_j->x, atom_i->rel_map,
                        atom_j->rel_map, nbr_pj->rel_box,
                        nbr_pj->dvec );

                if ( nbr_pj->d <= control->nonb_cut )
                {
                    flag = TRUE;

                    if ( nbr_pj->d <= control->nonb_sp_cut )
                    {
                        flag_sp = TRUE;
                    }
                }
            }

            if ( flag == TRUE )
            {
                type_j = system->atoms[j].type;
                sbp_j = &system->reax_param.sbp[type_j];
                twbp = &system->reax_param.tbp[type_i][type_j];
                r_ij = nbr_pj->d;

                val = Init_Charge_Matrix_Entry( system, control,
                            workspace, i, j, r_ij, OFF_DIAGONAL );
                val_flag = FALSE;

                for ( target = H->start[i]; target < Htop; ++target )
                {
                    if ( H->j[target] == j )
                    {
                        H->val[target] += val;
                        val_flag = TRUE;
                        break;
                    }
                }

                if ( val_flag == FALSE )
                {
                    H->j[Htop] = j;
                    H->val[Htop] = val;
                    ++Htop;
                }

                /* H_sp matrix entry */
                if ( flag_sp == TRUE )
                {
                    val_flag = FALSE;

                    for ( target = H_sp->start[i]; target < H_sp_top; ++target )
                    {
                        if ( H_sp->j[target] == j )
                        {
                            H_sp->val[target] += val;
                            val_flag = TRUE;
                            break;
                        }
                    }

                    if ( val_flag == FALSE )
                    {
                        H_sp->j[H_sp_top] = j;
                        H_sp->val[H_sp_top] = val;
                        ++H_sp_top;
                    }
                }

                /* hydrogen bond lists */
                if ( control->hbond_cut > 0.0 && (ihb == H_ATOM || ihb == H_BONDING_ATOM)
                        && nbr_pj->d <= control->hbond_cut )
                {
                    // fprintf( stderr, "%d %d\n", atom1, atom2 );
                    jhb = sbp_j->p_hbond;

                    if ( ihb == H_ATOM && jhb == H_BONDING_ATOM )
                    {
                        hbonds->hbond_list[ihb_top].nbr = j;
                        hbonds->hbond_list[ihb_top].scl = 1;
                        hbonds->hbond_list[ihb_top].ptr = nbr_pj;
                        ++ihb_top;
                        ++num_hbonds;
                    }
                    else if ( ihb == H_BONDING_ATOM && jhb == H_ATOM )
                    {
                        jhb_top = End_Index( workspace->hbond_index[j], hbonds );
                        hbonds->hbond_list[jhb_top].nbr = i;
                        hbonds->hbond_list[jhb_top].scl = -1;
                        hbonds->hbond_list[jhb_top].ptr = nbr_pj;
                        Set_End_Index( workspace->hbond_index[j], jhb_top + 1, hbonds );
                        ++num_hbonds;
                    }
                }

                /* uncorrected bond orders */
                if ( nbr_pj->d <= control->bond_cut )
                {
                    r2 = SQR( r_ij );

                    if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0 )
                    {
                        C12 = twbp->p_bo1 * POW( r_ij / twbp->r_s, twbp->p_bo2 );
                        BO_s = (1.0 + control->bo_cut) * EXP( C12 );
                    }
                    else
                    {
                        C12 = 0.0;
                        BO_s = 0.0;
                    }

                    if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 )
                    {
                        C34 = twbp->p_bo3 * POW( r_ij / twbp->r_p, twbp->p_bo4 );
                        BO_pi = EXP( C34 );
                    }
                    else
                    {
                        C34 = 0.0;
                        BO_pi = 0.0;
                    }

                    if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 )
                    {
                        C56 = twbp->p_bo5 * POW( r_ij / twbp->r_pp, twbp->p_bo6 );
                        BO_pi2 = EXP( C56 );
                    }
                    else
                    {
                        C56 = 0.0;
                        BO_pi2 = 0.0;
                    }

                    /* Initially BO values are the uncorrected ones, page 1 */
                    BO = BO_s + BO_pi + BO_pi2;

                    if ( BO >= control->bo_cut )
                    {
                        num_bonds += 2;
                        /****** bonds i-j and j-i ******/
                        ibond = &bonds->bond_list[btop_i];
                        btop_j = End_Index( j, bonds );
                        jbond = &bonds->bond_list[btop_j];

                        ibond->nbr = j;
                        jbond->nbr = i;
                        ibond->d = r_ij;
                        jbond->d = r_ij;
                        rvec_Copy( ibond->dvec, nbr_pj->dvec );
                        rvec_Scale( jbond->dvec, -1, nbr_pj->dvec );
                        ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
                        ivec_Scale( jbond->rel_box, -1, nbr_pj->rel_box );
                        ibond->dbond_index = btop_i;
                        jbond->dbond_index = btop_i;
                        ibond->sym_index = btop_j;
                        jbond->sym_index = btop_i;
                        ++btop_i;
                        Set_End_Index( j, btop_j + 1, bonds );

                        bo_ij = &ibond->bo_data;
                        bo_ij->BO = BO;
                        bo_ij->BO_s = BO_s;
                        bo_ij->BO_pi = BO_pi;
                        bo_ij->BO_pi2 = BO_pi2;
                        bo_ji = &jbond->bo_data;
                        bo_ji->BO = BO;
                        bo_ji->BO_s = BO_s;
                        bo_ji->BO_pi = BO_pi;
                        bo_ji->BO_pi2 = BO_pi2;

                        /* Bond Order page2-3, derivative of total bond order prime */
                        Cln_BOp_s = twbp->p_bo2 * C12 / r2;
                        Cln_BOp_pi = twbp->p_bo4 * C34 / r2;
                        Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2;

                        /* Only dln_BOp_xx wrt. dr_i is stored here, note that
                           dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
                        rvec_Scale( bo_ij->dln_BOp_s, -bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
                        rvec_Scale( bo_ij->dln_BOp_pi, -bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
                        rvec_Scale( bo_ij->dln_BOp_pi2,
                                -bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
                        rvec_Scale( bo_ji->dln_BOp_s, -1., bo_ij->dln_BOp_s );
                        rvec_Scale( bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
                        rvec_Scale( bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );

                        /* Only dBOp wrt. dr_i is stored here, note that
                           dBOp/dr_i = -dBOp/dr_j and all others are 0 */
                        rvec_Scale( bo_ij->dBOp, -(bo_ij->BO_s * Cln_BOp_s
                                    + bo_ij->BO_pi * Cln_BOp_pi
                                    + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
                        rvec_Scale( bo_ji->dBOp, -1., bo_ij->dBOp );

                        rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );