/*----------------------------------------------------------------------
  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 "reaxc_ffield.h"
#include <mpi.h>
#include <cctype>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "reaxc_defs.h"
#include "error.h"
#include "reaxc_tool_box.h"

char Read_Force_Field( FILE *fp, reax_interaction *reax,
                       control_params *control )
{
  char    *s;
  char   **tmp;
  char ****tor_flag;
  int      c, i, j, k, l, m, n, o, p, cnt;
  int lgflag = control->lgflag;
  int errorflag = 1;
  double     val;
  MPI_Comm comm;
  int me;
  comm = MPI_COMM_WORLD;
  MPI_Comm_rank(comm, &me);

  s = (char*) malloc(sizeof(char)*MAX_LINE);
  tmp = (char**) malloc(sizeof(char*)*MAX_TOKENS);
  for (i=0; i < MAX_TOKENS; i++)
    tmp[i] = (char*) malloc(sizeof(char)*MAX_TOKEN_LEN);

  /* reading first header comment */
  fgets( s, MAX_LINE, fp );

  /* line 2 is number of global parameters */
  fgets( s, MAX_LINE, fp );
  c = Tokenize( s, &tmp );

  /* reading the number of global parameters */
  n = atoi(tmp[0]);
  if (n < 1) {
    if (me == 0)
      control->error_ptr->warning( FLERR, "Number of globals in ffield file is 0. The file will not be read." );
    fclose(fp);
    free(s);
    free(tmp);
    return 1;
  }

  reax->gp.n_global = n;
  reax->gp.l = (double*) malloc(sizeof(double)*n);

  /* see reax_types.h for mapping between l[i] and the lambdas used in ff */
  for (i=0; i < n; i++) {
    fgets(s,MAX_LINE,fp);
    c = Tokenize(s,&tmp);

    val = (double) atof(tmp[0]);
    reax->gp.l[i] = val;
  }

  control->bo_cut    = 0.01 * reax->gp.l[29];
  control->nonb_low  = reax->gp.l[11];
  control->nonb_cut  = reax->gp.l[12];

  /* next line is number of atom types and some comments */
  fgets( s, MAX_LINE, fp );
  c = Tokenize( s, &tmp );
  reax->num_atom_types = atoi(tmp[0]);

  /* 3 lines of comments */
  fgets(s,MAX_LINE,fp);
  fgets(s,MAX_LINE,fp);
  fgets(s,MAX_LINE,fp);

  /* Allocating structures in reax_interaction */
  reax->sbp = (single_body_parameters*)
    scalloc(control->error_ptr,  reax->num_atom_types, sizeof(single_body_parameters), "sbp");
  reax->tbp = (two_body_parameters**)
    scalloc(control->error_ptr,  reax->num_atom_types, sizeof(two_body_parameters*), "tbp");
  reax->thbp= (three_body_header***)
    scalloc(control->error_ptr,  reax->num_atom_types, sizeof(three_body_header**), "thbp");
  reax->hbp = (hbond_parameters***)
    scalloc(control->error_ptr,  reax->num_atom_types, sizeof(hbond_parameters**), "hbp");
  reax->fbp = (four_body_header****)
    scalloc(control->error_ptr,  reax->num_atom_types, sizeof(four_body_header***), "fbp");
  tor_flag  = (char****)
    scalloc(control->error_ptr,  reax->num_atom_types, sizeof(char***), "tor_flag");

  for( i = 0; i < reax->num_atom_types; i++ ) {
    reax->tbp[i] = (two_body_parameters*)
      scalloc(control->error_ptr,  reax->num_atom_types, sizeof(two_body_parameters), "tbp[i]");
    reax->thbp[i]= (three_body_header**)
      scalloc(control->error_ptr,  reax->num_atom_types, sizeof(three_body_header*), "thbp[i]");
    reax->hbp[i] = (hbond_parameters**)
      scalloc(control->error_ptr,  reax->num_atom_types, sizeof(hbond_parameters*), "hbp[i]");
    reax->fbp[i] = (four_body_header***)
      scalloc(control->error_ptr,  reax->num_atom_types, sizeof(four_body_header**), "fbp[i]");
    tor_flag[i]  = (char***)
      scalloc(control->error_ptr,  reax->num_atom_types, sizeof(char**), "tor_flag[i]");

    for( j = 0; j < reax->num_atom_types; j++ ) {
      reax->thbp[i][j]= (three_body_header*)
        scalloc(control->error_ptr,  reax->num_atom_types, sizeof(three_body_header), "thbp[i,j]");
      reax->hbp[i][j] = (hbond_parameters*)
        scalloc(control->error_ptr,  reax->num_atom_types, sizeof(hbond_parameters), "hbp[i,j]");
      reax->fbp[i][j] = (four_body_header**)
        scalloc(control->error_ptr,  reax->num_atom_types, sizeof(four_body_header*), "fbp[i,j]");
      tor_flag[i][j]  = (char**)
        scalloc(control->error_ptr,  reax->num_atom_types, sizeof(char*), "tor_flag[i,j]");

      for (k=0; k < reax->num_atom_types; k++) {
        reax->fbp[i][j][k] = (four_body_header*)
          scalloc(control->error_ptr,  reax->num_atom_types, sizeof(four_body_header), "fbp[i,j,k]");
        tor_flag[i][j][k]  = (char*)
          scalloc(control->error_ptr,  reax->num_atom_types, sizeof(char), "tor_flag[i,j,k]");
      }
    }
  }

  reax->gp.vdw_type = 0;

  char errmsg[1024];

  for( i = 0; i < reax->num_atom_types; i++ ) {
    /* line one */
    fgets( s, MAX_LINE, fp );
    c = Tokenize( s, &tmp );

    /* Sanity checks */
    if (c == 2 && !lgflag)
        control->error_ptr->all(FLERR, "Force field file requires using 'lgvdw yes'");

    if (c < 9) {
      snprintf (errmsg, 1024, "Missing parameter(s) in line %s", s);
      control->error_ptr->all(FLERR, errmsg);
    }

    for( j = 0; j < (int)(strlen(tmp[0])); ++j )
      reax->sbp[i].name[j] = toupper( tmp[0][j] );

    val = atof(tmp[1]); reax->sbp[i].r_s        = val;
    val = atof(tmp[2]); reax->sbp[i].valency    = val;
    val = atof(tmp[3]); reax->sbp[i].mass       = val;
    val = atof(tmp[4]); reax->sbp[i].r_vdw      = val;
    val = atof(tmp[5]); reax->sbp[i].epsilon    = val;
    val = atof(tmp[6]); reax->sbp[i].gamma      = val;
    val = atof(tmp[7]); reax->sbp[i].r_pi       = val;
    val = atof(tmp[8]); reax->sbp[i].valency_e  = val;
    reax->sbp[i].nlp_opt = 0.5 * (reax->sbp[i].valency_e-reax->sbp[i].valency);

    /* line two */
    fgets( s, MAX_LINE, fp );
    c = Tokenize( s, &tmp );

    /* Sanity check */
    if (c < 8) {
      snprintf (errmsg, 1024, "Missing parameter(s) in line %s", s);
      control->error_ptr->all(FLERR, errmsg);
    }

    val = atof(tmp[0]); reax->sbp[i].alpha      = val;
    val = atof(tmp[1]); reax->sbp[i].gamma_w    = val;
    val = atof(tmp[2]); reax->sbp[i].valency_boc= val;
    val = atof(tmp[3]); reax->sbp[i].p_ovun5    = val;
    val = atof(tmp[4]);
    val = atof(tmp[5]); reax->sbp[i].chi        = val;
    val = atof(tmp[6]); reax->sbp[i].eta        = 2.0 * val;
    val = atof(tmp[7]); reax->sbp[i].p_hbond = (int) val;

    /* line 3 */
    fgets( s, MAX_LINE, fp );
    c = Tokenize( s, &tmp );

    /* Sanity check */
    if (c < 8) {
      snprintf (errmsg, 1024, "Missing parameter(s) in line %s", s);
      control->error_ptr->all(FLERR, errmsg);
    }

    val = atof(tmp[0]); reax->sbp[i].r_pi_pi    = val;
    val = atof(tmp[1]); reax->sbp[i].p_lp2      = val;
    val = atof(tmp[2]);
    val = atof(tmp[3]); reax->sbp[i].b_o_131    = val;
    val = atof(tmp[4]); reax->sbp[i].b_o_132    = val;
    val = atof(tmp[5]); reax->sbp[i].b_o_133    = val;
    val = atof(tmp[6]);
    val = atof(tmp[7]);

    /* line 4  */
    fgets( s, MAX_LINE, fp );
    c = Tokenize( s, &tmp );

    /* Sanity check */
    if (c < 8) {
      snprintf (errmsg, 1024, "Missing parameter(s) in line %s", s);
      control->error_ptr->all(FLERR, errmsg);
    }

    val = atof(tmp[0]); reax->sbp[i].p_ovun2    = val;
    val = atof(tmp[1]); reax->sbp[i].p_val3     = val;
    val = atof(tmp[2]);
    val = atof(tmp[3]); reax->sbp[i].valency_val= val;
    val = atof(tmp[4]); reax->sbp[i].p_val5     = val;
    val = atof(tmp[5]); reax->sbp[i].rcore2     = val;
    val = atof(tmp[6]); reax->sbp[i].ecore2     = val;
    val = atof(tmp[7]); reax->sbp[i].acore2     = val;

    /* line 5, only if lgvdw is yes */
    if (lgflag) {
      fgets( s, MAX_LINE, fp );
      c = Tokenize( s, &tmp );

      /* Sanity check */
      if (c > 2) {
        control->error_ptr->all(FLERR,"Force field file incompatible with 'lgvdw yes'");
      }

      val = atof(tmp[0]); reax->sbp[i].lgcij           = val;
      val = atof(tmp[1]); reax->sbp[i].lgre           = val;
    }

    if (reax->sbp[i].rcore2>0.01 && reax->sbp[i].acore2>0.01) { // Inner-wall
      if (reax->sbp[i].gamma_w>0.5) { // Shielding vdWaals
        if (reax->gp.vdw_type != 0 && reax->gp.vdw_type != 3) {
          if (errorflag && (me == 0)) {
            char errmsg[512];
            snprintf(errmsg, 512, "VdWaals-parameters for element %s "
                   "indicate inner wall+shielding, but earlier "
                   "atoms indicate different vdWaals-method. "
                   "This may cause division-by-zero errors. "
                   "Keeping vdWaals-setting for earlier atoms.",
                   reax->sbp[i].name);
            control->error_ptr->warning(FLERR,errmsg);
          }
          errorflag = 0;
        } else {
          reax->gp.vdw_type = 3;
        }
      } else {  // No shielding vdWaals parameters present
        if (reax->gp.vdw_type != 0 && reax->gp.vdw_type != 2) {
          if (me == 0) {
            char errmsg[512];
            snprintf(errmsg, 512, "VdWaals-parameters for element %s "
                   "indicate inner wall without shielding, but earlier "
                   "atoms indicate different vdWaals-method. "
                   "This may cause division-by-zero errors. "
                   "Keeping vdWaals-setting for earlier atoms.",
                   reax->sbp[i].name);
            control->error_ptr->warning(FLERR,errmsg);
          }
        } else {
          reax->gp.vdw_type = 2;
        }
      }
    } else { // No Inner wall parameters present
      if (reax->sbp[i].gamma_w>0.5) { // Shielding vdWaals
        if (reax->gp.vdw_type != 0 && reax->gp.vdw_type != 1) {
          if (me == 0) {
            char errmsg[512];
            snprintf(errmsg, 512, "VdWaals parameters for element %s "
                   "indicate shielding without inner wall, but earlier "
                   "elements indicate different vdWaals-method. "
                   "This may cause division-by-zero errors. "
                   "Keeping vdWaals-setting for earlier atoms.",
                   reax->sbp[i].name);
            control->error_ptr->warning(FLERR,errmsg);
          }
        } else {
          reax->gp.vdw_type = 1;
        }
      } else {
        char errmsg[256];
        snprintf(errmsg, 256, "Inconsistent vdWaals-parameters: "
                 "No shielding or inner-wall set for element %s",
                 reax->sbp[i].name);
        control->error_ptr->all(FLERR, errmsg);
      }
    }
  }

  /* Equate vval3 to valf for first-row elements (25/10/2004) */
  for( i = 0; i < reax->num_atom_types; i++ )
    if( reax->sbp[i].mass < 21 &&
        reax->sbp[i].valency_val != reax->sbp[i].valency_boc ) {
      if (me == 0) {
        char errmsg[256];
        snprintf(errmsg, 256, "Changed valency_val to valency_boc for %s",
               reax->sbp[i].name);
        control->error_ptr->warning(FLERR,errmsg);
      }
      reax->sbp[i].valency_val = reax->sbp[i].valency_boc;
    }

  /* next line is number of two body combination and some comments */
  fgets(s,MAX_LINE,fp);
  c=Tokenize(s,&tmp);

  if (c == 2 && !lgflag) {
      control->error_ptr->all(FLERR, "Force field file requires using 'lgvdw yes'");
    }

  l = atoi(tmp[0]);

  /* a line of comments */
  fgets(s,MAX_LINE,fp);

  for (i=0; i < l; i++) {
    /* line 1 */
    fgets(s,MAX_LINE,fp);
    c=Tokenize(s,&tmp);

    j = atoi(tmp[0]) - 1;
    k = atoi(tmp[1]) - 1;
    if ((c < 10) || (j < 0) || (k < 0))
      control->error_ptr->all(FLERR, "Inconsistent force field file");

    if (j < reax->num_atom_types && k < reax->num_atom_types) {

      val = atof(tmp[2]); reax->tbp[j][k].De_s      = val;
      reax->tbp[k][j].De_s      = val;
      val = atof(tmp[3]); reax->tbp[j][k].De_p      = val;
      reax->tbp[k][j].De_p      = val;
      val = atof(tmp[4]); reax->tbp[j][k].De_pp     = val;
      reax->tbp[k][j].De_pp     = val;
      val = atof(tmp[5]); reax->tbp[j][k].p_be1     = val;
      reax->tbp[k][j].p_be1     = val;
      val = atof(tmp[6]); reax->tbp[j][k].p_bo5     = val;
      reax->tbp[k][j].p_bo5     = val;
      val = atof(tmp[7]); reax->tbp[j][k].v13cor    = val;
      reax->tbp[k][j].v13cor    = val;

      val = atof(tmp[8]); reax->tbp[j][k].p_bo6     = val;
      reax->tbp[k][j].p_bo6     = val;
      val = atof(tmp[9]); reax->tbp[j][k].p_ovun1 = val;
      reax->tbp[k][j].p_ovun1 = val;

      /* line 2 */
      fgets(s,MAX_LINE,fp);
      c=Tokenize(s,&tmp);
      if (c < 8)
        control->error_ptr->all(FLERR, "Inconsistent force field file");

      val = atof(tmp[0]); reax->tbp[j][k].p_be2     = val;
      reax->tbp[k][j].p_be2     = val;
      val = atof(tmp[1]); reax->tbp[j][k].p_bo3     = val;
      reax->tbp[k][j].p_bo3     = val;
      val = atof(tmp[2]); reax->tbp[j][k].p_bo4     = val;
      reax->tbp[k][j].p_bo4     = val;
      val = atof(tmp[3]);

      val = atof(tmp[4]); reax->tbp[j][k].p_bo1     = val;
      reax->tbp[k][j].p_bo1     = val;
      val = atof(tmp[5]); reax->tbp[j][k].p_bo2     = val;
      reax->tbp[k][j].p_bo2     = val;
      val = atof(tmp[6]); reax->tbp[j][k].ovc       = val;
      reax->tbp[k][j].ovc       = val;

      val = atof(tmp[7]);
    }
  }

  for (i=0; i < reax->num_atom_types; i++)
    for (j=i; j < reax->num_atom_types; j++) {
      reax->tbp[i][j].r_s = 0.5 *
        (reax->sbp[i].r_s + reax->sbp[j].r_s);
      reax->tbp[j][i].r_s = 0.5 *
        (reax->sbp[j].r_s + reax->sbp[i].r_s);

      reax->tbp[i][j].r_p = 0.5 *
        (reax->sbp[i].r_pi + reax->sbp[j].r_pi);
      reax->tbp[j][i].r_p = 0.5 *
        (reax->sbp[j].r_pi + reax->sbp[i].r_pi);

      reax->tbp[i][j].r_pp = 0.5 *
        (reax->sbp[i].r_pi_pi + reax->sbp[j].r_pi_pi);
      reax->tbp[j][i].r_pp = 0.5 *
        (reax->sbp[j].r_pi_pi + reax->sbp[i].r_pi_pi);


      reax->tbp[i][j].p_boc3 =
        sqrt(reax->sbp[i].b_o_132 *
             reax->sbp[j].b_o_132);
      reax->tbp[j][i].p_boc3 =
        sqrt(reax->sbp[j].b_o_132 *
             reax->sbp[i].b_o_132);

      reax->tbp[i][j].p_boc4 =
        sqrt(reax->sbp[i].b_o_131 *
             reax->sbp[j].b_o_131);
      reax->tbp[j][i].p_boc4 =
        sqrt(reax->sbp[j].b_o_131 *
             reax->sbp[i].b_o_131);

      reax->tbp[i][j].p_boc5 =
        sqrt(reax->sbp[i].b_o_133 *
             reax->sbp[j].b_o_133);
      reax->tbp[j][i].p_boc5 =
        sqrt(reax->sbp[j].b_o_133 *
             reax->sbp[i].b_o_133);


      reax->tbp[i][j].D =
        sqrt(reax->sbp[i].epsilon *
             reax->sbp[j].epsilon);

      reax->tbp[j][i].D =
        sqrt(reax->sbp[j].epsilon *
             reax->sbp[i].epsilon);

      reax->tbp[i][j].alpha =
        sqrt(reax->sbp[i].alpha *
             reax->sbp[j].alpha);

      reax->tbp[j][i].alpha =
        sqrt(reax->sbp[j].alpha *
             reax->sbp[i].alpha);

      reax->tbp[i][j].r_vdW =
        2.0 * sqrt(reax->sbp[i].r_vdw * reax->sbp[j].r_vdw);

      reax->tbp[j][i].r_vdW =
        2.0 * sqrt(reax->sbp[j].r_vdw * reax->sbp[i].r_vdw);

      reax->tbp[i][j].gamma_w =
        sqrt(reax->sbp[i].gamma_w *
             reax->sbp[j].gamma_w);

      reax->tbp[j][i].gamma_w =
        sqrt(reax->sbp[j].gamma_w *
             reax->sbp[i].gamma_w);

      reax->tbp[i][j].gamma =
        pow(reax->sbp[i].gamma *
            reax->sbp[j].gamma,-1.5);

      reax->tbp[j][i].gamma =
        pow(reax->sbp[j].gamma *
            reax->sbp[i].gamma,-1.5);

      // additions for additional vdWaals interaction types - inner core

      reax->tbp[i][j].rcore = reax->tbp[j][i].rcore =
        sqrt( reax->sbp[i].rcore2 * reax->sbp[j].rcore2 );

      reax->tbp[i][j].ecore = reax->tbp[j][i].ecore =
        sqrt( reax->sbp[i].ecore2 * reax->sbp[j].ecore2 );

      reax->tbp[i][j].acore = reax->tbp[j][i].acore =
        sqrt( reax->sbp[i].acore2 * reax->sbp[j].acore2 );

      // additions for additional vdWalls interaction types lg correction

      reax->tbp[i][j].lgcij = reax->tbp[j][i].lgcij =
        sqrt( reax->sbp[i].lgcij * reax->sbp[j].lgcij );

      reax->tbp[i][j].lgre = reax->tbp[j][i].lgre = 2.0 * reax->gp.l[35] *
        sqrt( reax->sbp[i].lgre*reax->sbp[j].lgre );

    }

  fgets(s,MAX_LINE,fp);
  c=Tokenize(s,&tmp);
  l = atoi(tmp[0]);

  for (i=0; i < l; i++) {
    fgets(s,MAX_LINE,fp);
    c=Tokenize(s,&tmp);

    j = atoi(tmp[0]) - 1;
    k = atoi(tmp[1]) - 1;

    if ((c < (lgflag ? 9 : 8)) || (j < 0) || (k < 0))
      control->error_ptr->all(FLERR, "Inconsistent force field file");

    if (j < reax->num_atom_types && k < reax->num_atom_types)        {
      val = atof(tmp[2]);
      if (val > 0.0) {
        reax->tbp[j][k].D = val;
        reax->tbp[k][j].D = val;
      }

      val = atof(tmp[3]);
      if (val > 0.0) {
        reax->tbp[j][k].r_vdW = 2 * val;
        reax->tbp[k][j].r_vdW = 2 * val;
      }

      val = atof(tmp[4]);
      if (val > 0.0) {
        reax->tbp[j][k].alpha = val;
        reax->tbp[k][j].alpha = val;
      }

      val = atof(tmp[5]);
      if (val > 0.0) {
        reax->tbp[j][k].r_s = val;
        reax->tbp[k][j].r_s = val;
      }

      val = atof(tmp[6]);
      if (val > 0.0) {
        reax->tbp[j][k].r_p = val;
        reax->tbp[k][j].r_p = val;
      }

      val = atof(tmp[7]);
      if (val > 0.0) {
        reax->tbp[j][k].r_pp = val;
        reax->tbp[k][j].r_pp = val;
      }

      if (lgflag) {
        val = atof(tmp[8]);
        if (val >= 0.0) {
          reax->tbp[j][k].lgcij = val;
          reax->tbp[k][j].lgcij = val;
        }
      }
    }
  }

  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 )
        reax->thbp[i][j][k].cnt = 0;

  fgets( s, MAX_LINE, fp );
  c = Tokenize( s, &tmp );
  l = atoi( tmp[0] );

  for( i = 0; i < l; i++ ) {
    fgets(s,MAX_LINE,fp);
    c=Tokenize(s,&tmp);

    j = atoi(tmp[0]) - 1;
    k = atoi(tmp[1]) - 1;
    m = atoi(tmp[2]) - 1;
    if ((c < 10) || (j < 0) || (k < 0) || (m < 0))
      control->error_ptr->all(FLERR, "Inconsistent force field file");

    if (j < reax->num_atom_types && k < reax->num_atom_types &&
        m < reax->num_atom_types) {
      cnt = reax->thbp[j][k][m].cnt;
      reax->thbp[j][k][m].cnt++;
      reax->thbp[m][k][j].cnt++;

      val = atof(tmp[3]);
      reax->thbp[j][k][m].prm[cnt].theta_00 = val;
      reax->thbp[m][k][j].prm[cnt].theta_00 = val;

      val = atof(tmp[4]);
      reax->thbp[j][k][m].prm[cnt].p_val1 = val;
      reax->thbp[m][k][j].prm[cnt].p_val1 = val;

      val = atof(tmp[5]);
      reax->thbp[j][k][m].prm[cnt].p_val2 = val;
      reax->thbp[m][k][j].prm[cnt].p_val2 = val;

      val = atof(tmp[6]);
      reax->thbp[j][k][m].prm[cnt].p_coa1 = val;
      reax->thbp[m][k][j].prm[cnt].p_coa1 = val;

      val = atof(tmp[7]);
      reax->thbp[j][k][m].prm[cnt].p_val7 = val;
      reax->thbp[m][k][j].prm[cnt].p_val7 = val;

      val = atof(tmp[8]);
      reax->thbp[j][k][m].prm[cnt].p_pen1 = val;
      reax->thbp[m][k][j].prm[cnt].p_pen1 = val;

      val = atof(tmp[9]);
      reax->thbp[j][k][m].prm[cnt].p_val4 = val;
      reax->thbp[m][k][j].prm[cnt].p_val4 = val;
    }
  }

  /* clear all entries first */
  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 )
        for( m = 0; m < reax->num_atom_types; ++m ) {
          reax->fbp[i][j][k][m].cnt = 0;
          tor_flag[i][j][k][m] = 0;
        }

  /* next line is number of 4-body params and some comments */
  fgets( s, MAX_LINE, fp );
  c = Tokenize( s, &tmp );
  l = atoi( tmp[0] );

  for( i = 0; i < l; i++ ) {
    fgets( s, MAX_LINE, fp );
    c = Tokenize( s, &tmp );

    j = atoi(tmp[0]) - 1;
    k = atoi(tmp[1]) - 1;
    m = atoi(tmp[2]) - 1;
    n = atoi(tmp[3]) - 1;
    if ((c < 9) || (k < 0) || (m < 0))
      control->error_ptr->all(FLERR, "Inconsistent force field file");

    if (j >= 0 && n >= 0) { // this means the entry is not in compact form
      if (j < reax->num_atom_types && k < reax->num_atom_types &&
          m < reax->num_atom_types && n < reax->num_atom_types) {
        tor_flag[j][k][m][n] = 1;
        tor_flag[n][m][k][j] = 1;

        reax->fbp[j][k][m][n].cnt = 1;
        reax->fbp[n][m][k][j].cnt = 1;

        val = atof(tmp[4]);
        reax->fbp[j][k][m][n].prm[0].V1 = val;
        reax->fbp[n][m][k][j].prm[0].V1 = val;

        val = atof(tmp[5]);
        reax->fbp[j][k][m][n].prm[0].V2 = val;
        reax->fbp[n][m][k][j].prm[0].V2 = val;

        val = atof(tmp[6]);
        reax->fbp[j][k][m][n].prm[0].V3 = val;
        reax->fbp[n][m][k][j].prm[0].V3 = val;

        val = atof(tmp[7]);
        reax->fbp[j][k][m][n].prm[0].p_tor1 = val;
        reax->fbp[n][m][k][j].prm[0].p_tor1 = val;

        val = atof(tmp[8]);
        reax->fbp[j][k][m][n].prm[0].p_cot1 = val;
        reax->fbp[n][m][k][j].prm[0].p_cot1 = val;
      }
    } else { /* This means the entry is of the form 0-X-Y-0 */
      if (k < reax->num_atom_types && m < reax->num_atom_types)
        for( p = 0; p < reax->num_atom_types; p++ )
          for( o = 0; o < reax->num_atom_types; o++ ) {
            reax->fbp[p][k][m][o].cnt = 1;
            reax->fbp[o][m][k][p].cnt = 1;

            if (tor_flag[p][k][m][o] == 0) {
              reax->fbp[p][k][m][o].prm[0].V1 = atof(tmp[4]);
              reax->fbp[p][k][m][o].prm[0].V2 = atof(tmp[5]);
              reax->fbp[p][k][m][o].prm[0].V3 = atof(tmp[6]);
              reax->fbp[p][k][m][o].prm[0].p_tor1 = atof(tmp[7]);
              reax->fbp[p][k][m][o].prm[0].p_cot1 = atof(tmp[8]);
            }

            if (tor_flag[o][m][k][p] == 0) {
              reax->fbp[o][m][k][p].prm[0].V1 = atof(tmp[4]);
              reax->fbp[o][m][k][p].prm[0].V2 = atof(tmp[5]);
              reax->fbp[o][m][k][p].prm[0].V3 = atof(tmp[6]);
              reax->fbp[o][m][k][p].prm[0].p_tor1 = atof(tmp[7]);
              reax->fbp[o][m][k][p].prm[0].p_cot1 = atof(tmp[8]);
            }
          }
    }
  }

  /* next line is number of hydrogen bond params and some comments */
  fgets( s, MAX_LINE, fp );
  c = Tokenize( s, &tmp );
  l = atoi( tmp[0] );

  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 )
        reax->hbp[i][j][k].r0_hb = -1.0;

  for( i = 0; i < l; i++ ) {
    fgets( s, MAX_LINE, fp );
    c = Tokenize( s, &tmp );

    j = atoi(tmp[0]) - 1;
    k = atoi(tmp[1]) - 1;
    m = atoi(tmp[2]) - 1;
    if ((c < 7) || (j < 0) || (k < 0) || (m < 0))
      control->error_ptr->all(FLERR, "Inconsistent force field file");

    if (j < reax->num_atom_types && m < reax->num_atom_types) {
      val = atof(tmp[3]);
      reax->hbp[j][k][m].r0_hb = val;

      val = atof(tmp[4]);
      reax->hbp[j][k][m].p_hb1 = val;

      val = atof(tmp[5]);
      reax->hbp[j][k][m].p_hb2 = val;

      val = atof(tmp[6]);
      reax->hbp[j][k][m].p_hb3 = val;
    }
  }

  /* deallocate helper storage */
  for( i = 0; i < MAX_TOKENS; i++ )
    free( tmp[i] );
  free( tmp );
  free( s );


  /* deallocate tor_flag */
  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++ ) {
        free( tor_flag[i][j][k] );
      }
      free( tor_flag[i][j] );
    }
    free( tor_flag[i] );
  }
  free( tor_flag );

  // close file

  fclose(fp);

  return SUCCESS;
}
