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 "single_body_interactions.h"
#include "two_body_interactions.h"
#include "three_body_interactions.h"
#include "four_body_interactions.h"
#include "list.h"
#include "print_utils.h"
#include "system_props.h"
#include "charges.h"
#include "vector.h"


/* File scope variables */
static interaction_function Interaction_Functions[NO_OF_INTERACTIONS];


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


void Dummy_Interaction( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace,
        list **lists, output_controls *out_control )
{
}


void Init_Bonded_Force_Functions( control_params *control )
{
    Interaction_Functions[0] = Calculate_Bond_Orders;
    Interaction_Functions[1] = Bond_Energy;  //*/Dummy_Interaction;
    Interaction_Functions[2] = LonePair_OverUnder_Coordination_Energy;
    //*/Dummy_Interaction;
    Interaction_Functions[3] = Three_Body_Interactions; //*/Dummy_Interaction;
    Interaction_Functions[4] = Four_Body_Interactions;  //*/Dummy_Interaction;
    if ( control->hb_cut > 0.0 )
    {
        Interaction_Functions[5] = Hydrogen_Bonds; //*/Dummy_Interaction;
    }
    else
    {
        Interaction_Functions[5] = Dummy_Interaction;
    }
    Interaction_Functions[6] = Dummy_Interaction; //empty
    Interaction_Functions[7] = Dummy_Interaction; //empty
    Interaction_Functions[8] = Dummy_Interaction; //empty
    Interaction_Functions[9] = Dummy_Interaction; //empty
}


void Compute_Bonded_Forces( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace,
        list **lists, output_controls *out_control )
{

    int i;
    //real t_start, t_end, t_elapsed;

#ifdef 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

    /* Implement all the function calls as function pointers */
    for ( i = 0; i < NO_OF_INTERACTIONS; i++ )
    {
        (Interaction_Functions[i])( system, control, data, workspace,
                lists, out_control );

#if defined(DEBUG_FOCUS)
        fprintf( stderr, "f%d-", i );
#endif

#ifdef TEST_FORCES
        (Print_Interactions[i])(system, control, data, workspace,
                                lists, out_control);
#endif
    }
}


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

#ifdef 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, lists[FAR_NBRS], out_control );
    t_elapsed = Get_Timing_Info( t_start );
    data->timing.cm += t_elapsed;

#if defined(DEBUG_FOCUS)
    fprintf( stderr, "qeq - " );
#endif

    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(DEBUG_FOCUS)
    fprintf( stderr, "nonb forces - " );
#endif

#ifdef 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! */
void Compute_Total_Force( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace, list **lists )
{
    int i;
    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->select.bond_list[pj].nbr )
                {
                    if ( control->ensemble == NVE || control->ensemble == NVT || 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
        #pragma omp barrier

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


void Validate_Lists( static_storage *workspace, list **lists, int step, int n,
        int Hmax, int Htop, int num_bonds, int num_hbonds )
{
    int i, flag;
    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,
                     "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, "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->num_intrs - 2 )
    {
        workspace->realloc.bonds = 1;

        if ( End_Index(i, bonds) > bonds->num_intrs )
        {
            fprintf( stderr, "step%d-bondchk failed: i=%d end(i)=%d bond_end=%d\n",
                     step, flag, End_Index(i, bonds), bonds->num_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, "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->num_intrs - Start_Index(i, hbonds)) * DANGER_ZONE )
        {
            workspace->realloc.hbonds = 1;

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


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

    switch ( control->charge_method )
    {
    case QEQ_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->reaxprm.sbp[system->atoms[i].type].eta;
            break;
            default:
                fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
                exit( INVALID_INPUT );
            break;
        }
        break;

    case EE_CM:
        //TODO
        switch ( pos )
        {
            case OFF_DIAGONAL:
            break;
            case DIAGONAL:
            break;
            default:
                fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
                exit( INVALID_INPUT );
            break;
        }
        break;

    case ACKS2_CM:
        //TODO
        switch ( pos )
        {
            case OFF_DIAGONAL:
            break;
            case DIAGONAL:
            break;
            default:
                fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
                exit( INVALID_INPUT );
            break;
        }
        break;

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

    return ret;
}


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

    ret = 0.0;

    switch ( control->charge_method )
    {
    case QEQ_CM:
        switch ( pos )
        {
            case OFF_DIAGONAL:
                Tap = control->Tap7 * r_ij + control->Tap6;
                Tap = Tap * r_ij + control->Tap5;
                Tap = Tap * r_ij + control->Tap4;
                Tap = Tap * r_ij + control->Tap3;
                Tap = Tap * r_ij + control->Tap2;
                Tap = Tap * r_ij + control->Tap1;
                Tap = Tap * r_ij + control->Tap0;

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

                ret = ((i == j) ? 0.5 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
            break;

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

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

    case EE_CM:
        switch ( pos )
        {
            case OFF_DIAGONAL:
                if ( r_ij < control->r_cut && r_ij > 0.001 )
                {
                    Tap = control->Tap7 * r_ij + control->Tap6;
                    Tap = Tap * r_ij + control->Tap5;
                    Tap = Tap * r_ij + control->Tap4;
                    Tap = Tap * r_ij + control->Tap3;
                    Tap = Tap * r_ij + control->Tap2;
                    Tap = Tap * r_ij + control->Tap1;
                    Tap = Tap * r_ij + control->Tap0;

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

                    ret = Tap * EV_to_KCALpMOL / dr3gamij_3;
                }
            break;

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

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

    case ACKS2_CM:
        switch ( pos )
        {
            case OFF_DIAGONAL:
                if ( r_ij < control->r_cut && r_ij > 0.001 )
                {
                    Tap = control->Tap7 * r_ij + control->Tap6;
                    Tap = Tap * r_ij + control->Tap5;
                    Tap = Tap * r_ij + control->Tap4;
                    Tap = Tap * r_ij + control->Tap3;
                    Tap = Tap * r_ij + control->Tap2;
                    Tap = Tap * r_ij + control->Tap1;
                    Tap = Tap * r_ij + control->Tap0;

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

                    ret = Tap * EV_to_KCALpMOL / dr3gamij_3;
                }
            break;

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

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

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

    return ret;
}


static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
        control_params *control, list *far_nbrs,
        sparse_matrix * H, sparse_matrix * H_sp,
        int * Htop, int * H_sp_top )
{
    int i, j, pj;
    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:
            if ( (X_diag = (real*) calloc(system->N, sizeof(real))) == NULL )
            {
                fprintf( stderr, "not enough memory for charge matrix. terminating.\n" );
                exit( INSUFFICIENT_MEMORY );
            }

            H->start[system->N] = *Htop;
            H_sp->start[system->N] = *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;
            }

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

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

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

                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;

                for ( pj = Start_Index(i, far_nbrs); pj < End_Index(i, far_nbrs); ++pj )
                {
                    if ( far_nbrs->select.far_nbr_list[pj].d <= control->r_cut )
                    {
                        j = far_nbrs->select.far_nbr_list[pj].nbr;

                        xcut = ( system->reaxprm.sbp[ system->atoms[i].type ].b_s_acks2
                                + system->reaxprm.sbp[ system->atoms[j].type ].b_s_acks2 )
                            / 2.0;

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

                            H->j[*Htop] = system->N + j + 1;
                            H->val[*Htop] = MAX( 0.0, bond_softness );
                            *Htop = *Htop + 1;

                            H_sp->j[*H_sp_top] = system->N + j + 1;
                            H_sp->val[*H_sp_top] = MAX( 0.0, bond_softness );
                            *H_sp_top = *H_sp_top + 1;

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

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

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

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

            for ( i = system->N + 1; i < system->N_cm - 1; ++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 - 1];
                        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 - 1];
                        break;
                    }
                }
            }

            for ( i = system->N + 1; 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;

            free( X_diag );
            break;

        default:
            break;
    }
}


void Init_Forces( reax_system *system, control_params *control,
        simulation_data *data, static_storage *workspace,
        list **lists, output_controls *out_control )
{
    int i, j, pj;
    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;
    real r_ij, r2;
    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;
    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 = btop_j = 0;

    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->reaxprm.sbp[type_i]);
        ihb = ihb_top = -1;

        if ( control->hb_cut > 0 && (ihb = sbp_i->p_hbond) == 1 )
        {
            ihb_top = End_Index( workspace->hbond_index[i], hbonds );
        }

        for ( pj = start_i; pj < end_i; ++pj )
        {
            nbr_pj = &( far_nbrs->select.far_nbr_list[pj] );
            j = nbr_pj->nbr;
            atom_j = &(system->atoms[j]);

            flag = 0;
            flag_sp = 0;
            if ((data->step - data->prev_steps) % control->reneighbor == 0)
            {
                if ( nbr_pj->d <= control->r_cut )
                {
                    flag = 1;
                    if ( nbr_pj->d <= control->r_sp_cut )
                    {
                        flag_sp = 1;
                    }
                }
                else
                {
                    flag = 0;
                    flag_sp = 0;
                }
            }
            else if ((nbr_pj->d = Sq_Distance_on_T3(atom_i->x, atom_j->x, &(system->box),
                                                    nbr_pj->dvec)) <= SQR(control->r_cut))
            {
                if ( nbr_pj->d <= SQR(control->r_sp_cut))
                {
                    flag_sp = 1;
                }
                nbr_pj->d = SQRT( nbr_pj->d );
                flag = 1;
            }

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

                H->j[Htop] = j;
                H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, j, 
                        r_ij, OFF_DIAGONAL );
                ++Htop;

                /* H_sp matrix entry */
                if ( flag_sp )
                {
                    H_sp->j[H_sp_top] = j;
                    H_sp->val[H_sp_top] = H->val[Htop - 1];
                    ++H_sp_top;
                }

                /* hydrogen bond lists */
                if ( control->hb_cut > 0 && (ihb == 1 || ihb == 2) &&
                        nbr_pj->d <= control->hb_cut )
                {
                    // fprintf( stderr, "%d %d\n", atom1, atom2 );
                    jhb = sbp_j->p_hbond;
                    if ( ihb == 1 && jhb == 2 )
                    {
                        hbonds->select.hbond_list[ihb_top].nbr = j;
                        hbonds->select.hbond_list[ihb_top].scl = 1;
                        hbonds->select.hbond_list[ihb_top].ptr = nbr_pj;
                        ++ihb_top;
                        ++num_hbonds;
                    }
                    else if ( ihb == 2 && jhb == 1 )
                    {
                        jhb_top = End_Index( workspace->hbond_index[j], hbonds );
                        hbonds->select.hbond_list[jhb_top].nbr = i;
                        hbonds->select.hbond_list[jhb_top].scl = -1;
                        hbonds->select.hbond_list[jhb_top].ptr = nbr_pj;
                        Set_End_Index( workspace->hbond_index[j], jhb_top + 1, hbonds );
                        ++num_hbonds;
                    }
                }

                /* uncorrected bond orders */
                if ( far_nbrs->select.far_nbr_list[pj].d <= control->nbr_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
                    {
                        BO_s = 0.0;
                        C12 = 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
                    {
                        BO_pi = 0.0;
                        C34 = 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
                    {
                        BO_pi2 = 0.0;
                        C56 = 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->select.bond_list[btop_i] );
                        btop_j = End_Index( j, bonds );
                        jbond = &(bonds->select.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_ji = &( jbond->bo_data );
                        bo_ji->BO = bo_ij->BO = BO;
                        bo_ji->BO_s = bo_ij->BO_s = BO_s;
                        bo_ji->BO_pi = bo_ij->BO_pi = BO_pi;
                        bo_ji->BO_pi2 = bo_ij->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 );
                        rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp );

                        bo_ij->BO_s -= control->bo_cut;
                        bo_ij->BO -= control->bo_cut;
                        bo_ji->BO_s -= control->bo_cut;
                        bo_ji->BO -= control->bo_cut;
                        workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
                        workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
                        bo_ij->Cdbo = 0.0;
                        bo_ij->Cdbopi = 0.0;
                        bo_ij->Cdbopi2 = 0.0;
                        bo_ji->Cdbo = 0.0;
                        bo_ji->Cdbopi = 0.0;
                        bo_ji->Cdbopi2 = 0.0;

                        /*fprintf( stderr, "%d %d %g %g %g\n",
                          i+1, j+1, bo_ij->BO, bo_ij->BO_pi, bo_ij->BO_pi2 );*/

                        /*fprintf( stderr, "Cln_BOp_s: %f, pbo2: %f, C12:%f\n",
                          Cln_BOp_s, twbp->p_bo2, C12 );
                          fprintf( stderr, "Cln_BOp_pi: %f, pbo4: %f, C34:%f\n",
                          Cln_BOp_pi, twbp->p_bo4, C34 );
                          fprintf( stderr, "Cln_BOp_pi2: %f, pbo6: %f, C56:%f\n",
                          Cln_BOp_pi2, twbp->p_bo6, C56 );*/
                        /*fprintf(stderr, "pbo1: %f, pbo2:%f\n", twbp->p_bo1, twbp->p_bo2);
                          fprintf(stderr, "pbo3: %f, pbo4:%f\n", twbp->p_bo3, twbp->p_bo4);
                          fprintf(stderr, "pbo5: %f, pbo6:%f\n", twbp->p_bo5, twbp->p_bo6);
                          fprintf( stderr, "r_s: %f, r_p: %f, r_pp: %f\n",
                          twbp->r_s, twbp->r_p, twbp->r_pp );
                          fprintf( stderr, "C12: %g, C34:%g, C56:%g\n", C12, C34, C56 );*/

                        /*fprintf( stderr, "\tfactors: %g %g %g\n",
                          -(bo_ij->BO_s * Cln_BOp_s + bo_ij->BO_pi * Cln_BOp_pi +
                          bo_ij->BO_pi2 * Cln_BOp_pp),
                          -bo_ij->BO_pi * Cln_BOp_pi, -bo_ij->BO_pi2 * Cln_BOp_pi2 );*/
                        /*fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
                          bo_ij->dBOp[0], bo_ij->dBOp[1], bo_ij->dBOp[2] );
                          fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
                          bo_ij->dln_BOp_pi[0], bo_ij->dln_BOp_pi[1],
                          bo_ij->dln_BOp_pi[2] );
                          fprintf( stderr, "dBOpi2:\t[%g, %g, %g]\n\n",
                          bo_ij->dln_BOp_pi2[0], bo_ij->dln_BOp_pi2[1],
                          bo_ij->dln_BOp_pi2[2] );*/

                        Set_End_Index( j, btop_j + 1, bonds );
                    }
                }
            }
        }

        /* diagonal entry */
        H->j[Htop] = i;
        H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, i,
                r_ij, DIAGONAL );
        ++Htop;

        H_sp->j[H_sp_top] = i;
        H_sp->val[H_sp_top] = H->val[Htop - 1];
        ++H_sp_top;

        Set_End_Index( i, btop_i, bonds );
        if ( ihb == 1 )
        {
            Set_End_Index( workspace->hbond_index[i], ihb_top, hbonds );
        }
    }

    Init_Charge_Matrix_Remaining_Entries( system, control, far_nbrs,
            H, H_sp, &Htop, &H_sp_top );

    // mark the end of j list
    H->start[system->N_cm] = Htop;
    H_sp->start[system->N_cm] = H_sp_top;

    /* validate lists - decide if reallocation is required! */