text: /*----------------------------------------------------------------------
    PuReMD - Purdue ReaxFF Molecular Dynamics Program
  
    Copyright (2010) Purdue University
    Hasan Metin Aktulga, hmaktulga@lbl.gov
    Joseph Fogarty, jcfogart@mail.usf.edu
    Sagar Pandit, pandit@usf.edu
    Ananth Y Grama, ayg@cs.purdue.edu
  
    Please cite the related publication:
    H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama,
    "Parallel Reactive Molecular Dynamics: Numerical Methods and
    Algorithmic Techniques", Parallel Computing, in press.
  
    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 "pair_reaxc.h"
  #include "reaxc_hydrogen_bonds.h"
  #include "reaxc_bond_orders.h"
  #include "reaxc_list.h"
  #include "reaxc_valence_angles.h"
  #include "reaxc_vector.h"
  
  void Hydrogen_Bonds( reax_system *system, control_params *control,
                       simulation_data *data, storage *workspace,
                       reax_list **lists, output_controls * /*out_control*/ )
  {
    int  i, j, k, pi, pk;
    int  type_i, type_j, type_k;
    int  start_j, end_j, hb_start_j, hb_end_j;
    int  hblist[MAX_BONDS];
    int  itr, top;
    int  num_hb_intrs = 0;
    ivec rel_jk;
    double r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2;
    double e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3;
    rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk;
    rvec dvec_jk, force, ext_press;
    hbond_parameters *hbp;
    bond_order_data *bo_ij;
    bond_data *pbond_ij;
    far_neighbor_data *nbr_jk;
    reax_list *bonds, *hbonds;
    bond_data *bond_list;
    hbond_data *hbond_list;
  
    // tally variables
    double fi_tmp[3], fk_tmp[3], delij[3], delkj[3];
  
    bonds = (*lists) + BONDS;
    bond_list = bonds->select.bond_list;
    hbonds = (*lists) + HBONDS;
    hbond_list = hbonds->select.hbond_list;
  
    for( j = 0; j < system->n; ++j )
      if (system->reax_param.sbp[system->my_atoms[j].type].p_hbond == 1) {
        type_j     = system->my_atoms[j].type;
        start_j    = Start_Index(j, bonds);
        end_j      = End_Index(j, bonds);
        hb_start_j = Start_Index( system->my_atoms[j].Hindex, hbonds );
        hb_end_j   = End_Index( system->my_atoms[j].Hindex, hbonds );
        if (type_j < 0) continue;
  
        top = 0;
        for( pi = start_j; pi < end_j; ++pi )  {
          pbond_ij = &( bond_list[pi] );
          i = pbond_ij->nbr;
          type_i = system->my_atoms[i].type;
          if (type_i < 0) continue;
          bo_ij = &(pbond_ij->bo_data);
  
          if( system->reax_param.sbp[type_i].p_hbond == 2 &&
              bo_ij->BO >= HB_THRESHOLD )
            hblist[top++] = pi;
        }
  
        for( pk = hb_start_j; pk < hb_end_j; ++pk ) {
          /* set k's varibles */
          k = hbond_list[pk].nbr;
          type_k = system->my_atoms[k].type;
          if (type_k < 0) continue;
          nbr_jk = hbond_list[pk].ptr;
          r_jk = nbr_jk->d;
          rvec_Scale( dvec_jk, hbond_list[pk].scl, nbr_jk->dvec );
  
          for( itr = 0; itr < top; ++itr ) {
            pi = hblist[itr];
            pbond_ij = &( bonds->select.bond_list[pi] );
            i = pbond_ij->nbr;
  
            if (system->my_atoms[i].orig_id != system->my_atoms[k].orig_id) {
              bo_ij = &(pbond_ij->bo_data);
              type_i = system->my_atoms[i].type;
              if (type_i < 0) continue;
              hbp = &(system->reax_param.hbp[ type_i ][ type_j ][ type_k ]);
              if (hbp->r0_hb <= 0.0) continue;
              ++num_hb_intrs;
  
              Calculate_Theta( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
                               &theta, &cos_theta );
              /* the derivative of cos(theta) */
              Calculate_dCos_Theta( pbond_ij->dvec, pbond_ij->d, dvec_jk, r_jk,
                                    &dcos_theta_di, &dcos_theta_dj,
                                    &dcos_theta_dk );
  
              /* hyrogen bond energy*/
              sin_theta2 = sin( theta/2.0 );
              sin_xhz4 = SQR(sin_theta2);
              sin_xhz4 *= sin_xhz4;
              cos_xhz1 = ( 1.0 - cos_theta );
              exp_hb2 = exp( -hbp->p_hb2 * bo_ij->BO );
              exp_hb3 = exp( -hbp->p_hb3 * ( hbp->r0_hb / r_jk +
                                             r_jk / hbp->r0_hb - 2.0 ) );
  
              data->my_en.e_hb += e_hb =
                hbp->p_hb1 * (1.0 - exp_hb2) * exp_hb3 * sin_xhz4;
  
              CEhb1 = hbp->p_hb1 * hbp->p_hb2 * exp_hb2 * exp_hb3 * sin_xhz4;
              CEhb2 = -hbp->p_hb1/2.0 * (1.0 - exp_hb2) * exp_hb3 * cos_xhz1;
              CEhb3 = -hbp->p_hb3 *
                (-hbp->r0_hb / SQR(r_jk) + 1.0 / hbp->r0_hb) * e_hb;
  
              /* hydrogen bond forces */
              bo_ij->Cdbo += CEhb1; // dbo term
  
              if (control->virial == 0) {
                // dcos terms
                rvec_ScaledAdd( workspace->f[i], +CEhb2, dcos_theta_di );
                rvec_ScaledAdd( workspace->f[j], +CEhb2, dcos_theta_dj );
                rvec_ScaledAdd( workspace->f[k], +CEhb2, dcos_theta_dk );
                // dr terms
                rvec_ScaledAdd( workspace->f[j], -CEhb3/r_jk, dvec_jk );
                rvec_ScaledAdd( workspace->f[k], +CEhb3/r_jk, dvec_jk );
              }
              else {
                rvec_Scale( force, +CEhb2, dcos_theta_di ); // dcos terms
                rvec_Add( workspace->f[i], force );
                rvec_iMultiply( ext_press, pbond_ij->rel_box, force );
                rvec_ScaledAdd( data->my_ext_press, 1.0, ext_press );
  
                rvec_ScaledAdd( workspace->f[j], +CEhb2, dcos_theta_dj );
  
                ivec_Scale( rel_jk, hbond_list[pk].scl, nbr_jk->rel_box );
                rvec_Scale( force, +CEhb2, dcos_theta_dk );
                rvec_Add( workspace->f[k], force );
                rvec_iMultiply( ext_press, rel_jk, force );
                rvec_ScaledAdd( data->my_ext_press, 1.0, ext_press );
                // dr terms
                rvec_ScaledAdd( workspace->f[j], -CEhb3/r_jk, dvec_jk );
  
                rvec_Scale( force, CEhb3/r_jk, dvec_jk );
                rvec_Add( workspace->f[k], force );
                rvec_iMultiply( ext_press, rel_jk, force );
                rvec_ScaledAdd( data->my_ext_press, 1.0, ext_press );
              }
  
              /* tally into per-atom virials */
              if (system->pair_ptr->vflag_atom || system->pair_ptr->evflag) {
                rvec_ScaledSum( delij, 1., system->my_atoms[j].x,
                                      -1., system->my_atoms[i].x );
                rvec_ScaledSum( delkj, 1., system->my_atoms[j].x,
                                       -1., system->my_atoms[k].x );
  
                rvec_Scale(fi_tmp, CEhb2, dcos_theta_di);
                rvec_Scale(fk_tmp, CEhb2, dcos_theta_dk);
                rvec_ScaledAdd(fk_tmp, CEhb3/r_jk, dvec_jk);
  
                system->pair_ptr->ev_tally3(i,j,k,e_hb,0.0,fi_tmp,fk_tmp,delij,delkj);
              }
            }
          }
        }
      }
  }
  