/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   Contributing author: Matthew Wander
------------------------------------------------------------------------- */

#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_vmm.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

#define MAXLINE 1024
#define DELTA 4

/* ---------------------------------------------------------------------- */

PairVMM::PairVMM(LAMMPS *lmp) : Pair(lmp)
{
  single_enable = 0;
  restartinfo = 0;
  one_coeff = 1;
  manybody_flag = 1;
  
  nelements = 0;
  elements = NULL;
  nparams = maxparam = 0;
  params = NULL;
  elem2param = NULL;
//  bonds= NULL;
  ewaldflag = pppmflag = 1;
}

/* ----------------------------------------------------------------------
   check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */

PairVMM::~PairVMM()
{
  if (copymode) return;

  if (elements)
    for (int i = 0; i < nelements; i++) delete [] elements[i];
  delete [] elements;
  memory->destroy(params);
  memory->destroy(elem2param);

  if (allocated) {
    memory->destroy(setflag);
//    memory->destroy(cutsq);
    delete [] map;
  }
}

/* ---------------------------------------------------------------------- */

void PairVMM::compute(int eflag, int vflag)
{
  int i,j,k,l,ii,jj,kk,nn,inum,jnum,ncoordi,ncoordj;
  int itype,jtype,ktype,ltype,ijparam,jiparam,jkparam,kjparam,idebug,is13bonded;
  tagint itag,jtag, ktag, ltag;
  double xtmp,ytmp,ztmp,delx,dely,delz,evdwl;
  double forceix,forceiy,forceiz,forcejx,forcejy,forcejz;
//  double ebond, evec;
  double rsq;
  int *ilist,*jlist,*numneigh,**firstneigh;
  
// Here are the new VMM variables. 
  double r, rjk, demu1,fmu1,bfd,f2d,stotal,stotalcorr;
  double sij,sjk,skl,jtemp;
  double vvsum,vvsideal,vvx,vvy,vvz,kvvs,vvdel;
  int mtype;

  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

  double **x = atom->x;
  double **f = atom->f;
  tagint *tag = atom->tag;
  int *type = atom->type;
  int nlocal = atom->nlocal;
  int newton_pair = force->newton_pair;

  inum = list->inum;
  ilist = list->ilist;
  numneigh = list->numneigh;
  firstneigh = list->firstneigh;

 // loop over full neighbor list of my atoms
  jnum=0;
  is13bonded=0;
  stotal=0;
  stotalcorr=0;

// set idebug here
//  idebug=1 debug mode
//  idebug=0 no debugging comments
    idebug=1;
///////////////////////////////////////////////////////
   
  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    itag = tag[i];
    itype = map[type[i]];
    xtmp = x[i][0];
    ytmp = x[i][1];
    ztmp = x[i][2];  

   /// zero the energy here. Yes this is different than otherwise but it works correctly and makes sense. 
    evdwl = 0.0;

    jlist = firstneigh[i];
    jnum = numneigh[i];
    
    if (idebug==1) {
      printf("------------------------------------------------------ \n");  
      printf("ii = %i \n",ii);
      printf("itype  = %i \n",itype);
      printf("jnum = %i \n",jnum);
    }
    
    // We need to hold informaiton on every atom bound to our ith atom. Here is where those variables are declared.
    double ssort[jnum],stotalj[jnum],stotaljcorr[jnum],vvxj[jnum],vvyj[jnum],vvzj[jnum], // sj[jnum],
                 kvvsj[jnum],vvsidealj[jnum],vvsumj[jnum],delxj[jnum],delyj[jnum],delzj[jnum],
                 sideali[jnum],sidealj[jnum],ssortj[jnum],saltconfig[jnum],sjjcorr[jj];  
                 double sjj[jnum][jnum],sjjexcess[jnum][jnum];

    int jjthelement[jnum][jnum], jthelement[jnum],kthelement[jnum];
        
    // zero zero zero and zero again. damn seg faults.  
    demu1=0;
    fmu1=0;
    bfd=0;
    stotal=0;
    
    jtemp=0;
    vvx=0;
    vvy=0;
    vvz=0;
    ncoordi=0;
    vvsum=0;
    vvdel=0;
    
    for (jj = 0; jj < jnum; jj++) {
      // sj[jj]=0;
      ssort[jj]=0;
      stotalj[jj]=0;
      stotaljcorr[jj]=0;
      //sjexcess[jj]=0;
      sidealj[jj]=0;
      vvxj[jj]=0;
      vvyj[jj]=0;
      vvzj[jj]=0;
      kvvsj[jj]=0;
      vvsidealj[jj]=0;
      jthelement[jj]=0;
      kthelement[jj]=0;
      delxj[jj]=0;
      delyj[jj]=0;
      delzj[jj]=0;
      sideali[jj]=0;
      sidealj[jj]=0;
      ssortj[jj]=0;
      sjjcorr[jj]=0;
      
      for (kk = 0; kk < jnum; kk++) {
          sjj[jj][kk]=0;
          // stotalk[jj][kk]=0;
          sjjexcess[jj][kk]=0;
          jjthelement[jj][kk]=0;
      } 
    }
    
    for (jj = 0; jj < jnum; jj++) {

      j = jlist[jj];
      j &= NEIGHMASK;
      jtag = tag[j];
      

      // zero some more
      sij=0;    
      
      jtype = map[type[j]];      

      ijparam = elem2param[itype][jtype];
      jiparam = elem2param[jtype][itype];

      jthelement[jj]=params[ijparam].jelement;
      jjthelement[jj][jj]=params[jiparam].jelement;

      if (params[ijparam].force_switch >= 0) {
        delx = xtmp - x[j][0];
        dely = ytmp - x[j][1];
        delz = ztmp - x[j][2];

        if (idebug==1) {
          printf("jj = %i \n",jj);
          printf("jtype  = %i \n",jtype);
        }    
           
        rsq = delx*delx + dely*dely + delz*delz;
        r = sqrt(rsq);
        
        // here is where we calcualte the bond valence for each bond. We also compute an "average" b value for that bond.
        if (r < params[ijparam].cut) {
          sij = bondvalence(&params[ijparam], r);

        }
        // sj[jj] = sij;
        stotal +=sij;
        stotalj[jj]=sij;
        sjj[jj][jj]=sij;
      }     
      
      //accumulate the paramters of our energy expression.

      if (idebug==1) {
        printf("Before loop k \n"); 
        printf("itype  = %i \n",itype);
        printf("jtype  = %i \n",jtype);
        printf("sij = %f \n",sij);
        //printf("sj[jj] = %f \n",sj[jj]);
        printf("stotal = %f \n",stotal);
        //printf("stotalj[jj] = %f \n",stotalj[jj]);
        //printf("ssort[jj] = %f \n",ssort[jj]);
        //printf("ssortj[jj] = %f \n",ssortj[jj]);
        printf("sjj[jj][jj] = %f \n",sjj[jj][jj]);
        
      }
      
      for (kk = 0; kk < jnum; kk++) {
        if (kk==jj) continue;  
        
        k = jlist[kk];
        k &= NEIGHMASK;
//        ktag = tag[k];
   
      // zero some more
        sjk=0;
        kthelement[kk]=0;      
        
        ktype = map[type[k]];      
        jkparam = elem2param[jtype][ktype];
        kjparam = elem2param[ktype][jtype];
        
        jjthelement[jj][kk]=params[jkparam].jelement;

        if (params[jkparam].force_switch >= 0) {
          delx = x[j][0] - x[k][0];
          dely = x[j][1] - x[k][1];
          delz = x[j][2] - x[k][2];
    
          rsq = delx*delx + dely*dely + delz*delz;
          r = sqrt(rsq);

//        printf("jkparam = %i \n",jkparam);
//        printf("r = %f \n",r);
          if (r < params[jkparam].cut) {
      // here is where we calcualte the bond valence for each bond. We also compute an "average" b value for that bond.
            sjk = bondvalence(&params[jkparam], r);
          }
            stotalj[jj]+=sjk;
            sjj[jj][kk] = sjk;
            // stotalk[jj][kk]=sjk;
          
        }

        if (idebug==1) {
          printf ("Inside loop k. \n");
          printf("jtype  = %i \n",jtype);
          printf("ktype  = %i \n",ktype);
          printf("sjk = %f \n",sjk);
          printf("sjj[jj][kk] = %f \n",sjj[jj][kk]);
        }
      }

      if (idebug==1) {
          printf("stotalj[jj] = %f \n",stotalj[jj]);
      }

      // for (kk = 0; kk < jnum; kk++) {
      //     ssortj[kk]=sjj[jj][kk];
      //     if (idebug==1) {
      //       printf ("Finish loop k. \n");
      //       printf("ssortj[kk] = %f \n",ssortj[kk]);
      //     }
      // }
      
      if (idebug==1) {
          printf ("Inbetween loops \n");
//    printf("vvsidealj[jj] = %f \n", vvsidealj[jj]);
//    printf("kvvsj[jj] = %f \n", kvvsj[jj]);
//    printf("ithelement = %i \n", ithelement);
//    printf("jthelement = %i \n", jthelement[jj]);
//    printf("ssort[0] = %f \n", ssort[0]);
//    printf("ssort[1] = %f \n", ssort[1]);
//    printf("start second loop here. \n");
//    printf("vvsideal = %f \n", vvsideal);
//    printf("kvvs = %f \n", kvvs);
      }
    }  

/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */    
// Start second (overbonding correction) loop here
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
    stotalcorr=stotal;
    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj];
      j &= NEIGHMASK;
      
      jtype = map[type[j]];   
      ijparam = elem2param[itype][jtype];
      jiparam = elem2param[jtype][itype]; 

      if ((stotal > params[ijparam].ideal_valence) && (stotalj[jj] > params[jiparam].ideal_valence) && (sjj[jj][jj]>0)) {
        sjjexcess[jj][jj]=std::min(std::min((stotalj[jj] - params[jiparam].ideal_valence), (stotal - params[ijparam].ideal_valence)),sjj[jj][jj]);
        //sjjexcess[jj][jj]=sjexcess[jj];

        if (idebug==1) {
            printf("In correction loop j\n");
            printf("itype  = %i \n",itype);
            printf("jtype  = %i \n",jtype);
            printf("stotal = %f \n",stotalj[jj]);
            printf("stotalj[jj] = %f \n",stotalj[jj]);
            printf("i ideal valence = %f \n",params[ijparam].ideal_valence);
            printf("j ideal valence = %f \n",params[jiparam].ideal_valence);
            printf("sjjexcess[jj][jj] = %f \n",sjjexcess[jj][jj]);
        } 
      }
    
      for (kk = 0; kk < jnum; kk++) {
        if (kk==jj) continue;
        
        k = jlist[kk];
        k &= NEIGHMASK;
        ktype = map[type[k]];      

        jkparam = elem2param[jtype][ktype];
        kjparam = elem2param[ktype][jtype];

        if ((stotalj[jj] > params[jkparam].ideal_valence) && (stotalj[kk] > params[kjparam].ideal_valence) && (sjj[jj][kk]>0)) {
          sjjexcess[jj][kk]=std::min(std::min((stotalj[jj] - params[jkparam].ideal_valence), (stotalj[kk] - params[kjparam].ideal_valence)),sjj[jj][kk]); 

          if (idebug==1) {
            printf("In correction loop k\n");
            printf("itype  = %i \n",jtype);
            printf("jtype  = %i \n",ktype);
            printf("stotalj[jj] = %f \n",stotalj[jj]);
            printf("stotalj[kk] = %f \n",stotalj[kk]);
            printf("sjj[jj][kk] = %f \n",sjj[jj][kk]);
            printf("sjjexcess[jj][kk] = %f \n",sjjexcess[jj][kk]);
          }          
        }
      }
      stotalcorr -= sjjexcess[jj][jj];
      ssort[jj] = sjj[jj][jj]-sjjexcess[jj][jj]; 
      stotaljcorr[jj] = stotalj[jj]; 
      
      for (kk = 0; kk < jnum; kk++) {
        stotaljcorr[jj] -= sjjexcess[jj][kk];
      }
      
      if (idebug==1) {
          printf("After correction loop k\n");
          printf("stotalj[jj] = %f \n",stotalj[jj]);
          printf("stotaljcorr[jj] = %f \n",stotaljcorr[jj]);
          printf("ssort[jj] = %f \n",ssort[jj]);
      }
    }
    if (idebug==1) {
      printf("After correction loop j \n");
      printf("stotalcorr = %f \n",stotalcorr);
    }



  
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */    
// Start third (energy) loop here
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */

    f[i][0] = 0.0; 
    f[i][1] = 0.0; 
    f[i][2] = 0.0;


    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj];
      j &= NEIGHMASK;
      // jtag = tag[j];

      jtype = map[type[j]];       

      delx = xtmp - x[j][0];
      dely = ytmp - x[j][1];
      delz = ztmp - x[j][2];
      r = sqrt(delx*delx + dely*dely + delz*delz);

      ijparam = elem2param[itype][jtype];
      jiparam = elem2param[jtype][itype];
      forceix=0;
      forceiy=0;
      forceiz=0;
      forcejx=0;
      forcejy=0;
      forcejz=0;


      if (idebug==1) {
        printf("Starting second J loop. \n");
        printf("jj = %i \n",jj);
        printf("jtype  = %i \n",jtype);
        printf("evdw = %f \n",evdwl);
        printf("stotalcorr = %f \n",stotalcorr);
        printf("stotaljcorr[jj] = %f \n",stotaljcorr[jj]);
        printf("sjj[jj][jj] = %f \n",sjj[jj][jj]);

//      printf("vvsum = %f \n",vvsum);
//      printf("vvsideal = %f \n",vvsideal);
//      printf("kvvs = %f \n",kvvs);
//      printf("vvx = %f \n",vvx);
//      printf("vvy = %f \n",vvy);
//     printf("vvz = %f \n",vvz);

      }      
 


      for (kk = 0; kk < jnum; kk++) {
        ssortj[kk]=sjj[jj][kk]-sjjexcess[jj][kk];
        sjjcorr[kk]=sjj[jj][kk]-sjjexcess[jj][kk];
        kthelement[kk]=jjthelement[jj][kk];
      }

      if (idebug==1) {
        printf("Starting second K loop. \n");
        printf("kk = %i \n",kk);
        printf("ktype  = %i \n",ktype);
      }
              
       // Call the coordination number routines of both i j
      ncoordi = ncoord(jnum, ssort, stotalcorr, jthelement);  //note jth and kth elements never assigned.
      ncoordj = ncoord(jnum, ssortj, stotaljcorr[jj], kthelement);
      
      if (idebug==1) {
        printf("ncoordi = %i \n",ncoordi);
        printf("ncoordj = %i \n",ncoordj);
      } 
      
      // Generate generic ideal configuration based on coordination number. 
      idealconfiguration(ncoordi, jnum, params[ijparam].ideal_valence, ssort, stotalcorr, params[ijparam].ielflag, 
                         jthelement, sideali);
      idealconfiguration(ncoordj, jnum, params[jiparam].ideal_valence, ssortj, stotaljcorr[jj], params[jiparam].ielflag, 
                         kthelement, sidealj);

      if (idebug==1){
        for (kk = 0; kk < jnum; kk++) {
            printf("sideali[ %i ] = %f \n",kk,sideali[kk]);  
            printf("sidealj[ %i ] = %f \n",kk,sidealj[kk]);          
        } 
      }
      
      vectorvalence(&params[jiparam], jnum, ssortj, sjjcorr, jj, jlist, ncoordj, xtmp, ytmp,ztmp, kthelement, stotalj[jj], vvsumj[jj], 
                  vvsidealj[jj], kvvsj[jj]);  // fix this. 

      /// This is where we compute the 13bonbed identifier.

      is13bonded=0;

      for (kk = 0; kk < jnum; kk++) {
        if (kk!=jj) {
        // note sjj[kk,kk] in this notation is actually sik.
          if (((sjj[kk][kk]-sjjexcess[kk][kk]) > 1/double(ncoordi*2)) && ((sjj[jj][kk]-sjjexcess[jj][kk]) > 1/double(ncoordj*2))) {             ///fix this for ideal bonds. Is that even possible?
            is13bonded=1;
          }
        }
      }
      
      if (idebug==1){
        printf("is13bonded = %i \n",is13bonded);          
      }
      
      // Now add up the energy for each atom and compute the forces.
      if (params[ijparam].force_switch == -10) {

        nonbondedenergy(&params[ijparam],r,eflag,is13bonded,evdwl,delx,dely,delz,forceix,forceiy,forceiz,idebug);

        f[i][0] += forceix;
        f[i][1] += forceiy;
        f[i][2] += forceiz;
        if (idebug==1) {
          printf("non bonded energy is = %f \n",evdwl);
        }
      }
      else {
        valencemonopole(&params[ijparam],params[jiparam].ideal_valence,r,jnum,eflag,is13bonded,evdwl,stotalcorr,stotaljcorr[jj],sjj[jj][jj]-sjjexcess[jj][jj],
              ssort,ssortj,sideali,sidealj,ncoordi,ncoordj,delxj,delyj,delzj,delx,dely,delz,forceix,forceiy,forceiz,i,j,idebug);
      
        f[i][0] += forceix;
        f[i][1] += forceiy;
        f[i][2] += forceiz;
        if (idebug==1) {
          printf("bonded energy is = %f \n",evdwl);
        }
      }
      
      // Now add up the energy for each atom and compute the forces.
      // Not sj is the ith jth bond so it is not part of the other matrix.

      if (idebug ==1)
      {    
        printf ("after two body before three body\n");  
        printf("stotalj[jj] = %f \n",stotalj[jj]);
        printf("fx-i = %f \n",f[i][0]);
        printf("fy-i = %f \n",f[i][1]);
        printf("fz-i = %f \n",f[i][2]);
        printf("vvsumj[jj] = %f \n",vvsumj[jj]);
        printf("vvsidealj[jj] = %f \n",vvsidealj[jj]);
        printf("kvvsj[jj] = %f \n",kvvsj[jj]);
        printf("vvxj[jj] = %f \n",vvxj[jj]);
        printf("vvyj[jj] = %f \n",vvyj[jj]);
        printf("vvzj[jj] = %f \n",vvzj[jj]);
        printf("evdwl = %f \n",evdwl); 
      }


      if (kvvsj[jj] > 0) {                
        valencedipole(&params[jiparam],r,eflag,evdwl,stotalj[jj],sjj[jj][jj]-sjjexcess[jj][jj],vvsumj[jj],vvsidealj[jj],vvxj[jj],vvyj[jj],vvzj[jj],kvvsj[jj],
              -delx,-dely,-delz,forcejx,forcejy,forcejz,idebug);

        f[i][0] -= forcejx;
        f[i][1] -= forcejy;
        f[i][2] -= forcejz;
        if (idebug==1) {
          printf("valence dipole energy is = %f \n",evdwl);
        }
      }
    
      if (idebug ==1) {
        printf ("after three body\n");  
        printf("fx-i = %f \n",f[i][0]);
        printf("fy-i = %f \n",f[i][1]);
        printf("fz-i = %f \n",f[i][2]);
//      printf("forceix = %f \n",forceix); 
//      printf("forceiy = %f \n",forceiy); 
//      printf("forceiz = %f \n",forceiz); 
//      printf("nlocal = %i \n",nlocal); 
//      printf("newton_pair = %i \n",newton_pair); 
        printf("final energy = %f \n",evdwl); 
//      printf("forcejx = %f \n",forcejx); 
//      printf("forcejy = %f \n",forcejy); 
//      printf("forcejz = %f \n",forcejz); 
      }

    }
     
    if (evflag) ev_tally_xyz_full(i, evdwl, 0.0, f[i][0], f[i][1], f[i][2], 1.0, 1.0, 1.0);
  }
//  if (vflag_fdotr) virial_fdotr_compute();
}

/* ---------------------------------------------------------------------- */

void PairVMM::allocate()
{
  allocated = 1;
  int n = atom->ntypes;

  memory->create(setflag,n+1,n+1,"pair:setflag");
  memory->create(cutsq,n+1,n+1,"pair:cutsq"); // While this line is unecessary removing it causes the program to segfault

  map = new int[n+1];
}

/* ----------------------------------------------------------------------
   global settings
------------------------------------------------------------------------- */

void PairVMM::settings(int narg, char **arg)
{
  if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}

/* ----------------------------------------------------------------------
   set coeffs for one or more type pairs
------------------------------------------------------------------------- */

void PairVMM::coeff(int narg, char **arg)
{
  int i,j,n;

  if (!allocated) allocate();

  if (narg != 3 + atom->ntypes)
    error->all(FLERR,"Incorrect args for pair coefficients");

  // insure I,J args are * *

  if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
    error->all(FLERR,"Incorrect args for pair coefficients");

  // read args that map atom types to elements in potential file
  // map[i] = which element the Ith atom type is, -1 if NULL
  // nelements = # of unique elements
  // elements = list of element names

  if (elements) {
    for (i = 0; i < nelements; i++) delete [] elements[i];
    delete [] elements;
  }
  elements = new char*[atom->ntypes];
  for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;

  nelements = 0;
  for (i = 3; i < narg; i++) {
    if (strcmp(arg[i],"NULL") == 0) {
      map[i-2] = -1;
      continue;
    }
    for (j = 0; j < nelements; j++)
      if (strcmp(arg[i],elements[j]) == 0) break;
    map[i-2] = j;
    if (j == nelements) {
      n = strlen(arg[i]) + 1;
      elements[j] = new char[n];
      strcpy(elements[j],arg[i]);
      nelements++;
    }
  }

  // read potential file and initialize potential parameters

  read_file(arg[2]);
  setup();

  // clear setflag since coeff() called once with I,J = * *

  n = atom->ntypes;
  for (int i = 1; i <= n; i++)
    for (int j = i; j <= n; j++)
      setflag[i][j] = 0;

  // set setflag i,j for type pairs where both are mapped to elements

  int count = 0;
  for (int i = 1; i <= n; i++)
    for (int j = i; j <= n; j++)
      if (map[i] >= 0 && map[j] >= 0) {
        setflag[i][j] = 1;
        count++;
      }

  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}

/* ----------------------------------------------------------------------
   init specific to this pair style
------------------------------------------------------------------------- */

void PairVMM::init_style()
{
  if (atom->tag_enable == 0)
    error->all(FLERR,"Pair style Valence Multipole Model requires atom IDs");
  if (force->newton_pair == 0)
    error->all(FLERR,"Pair style Valence Multipole Model requires newton pair on");

  // need a full neighbor list

  int irequest = neighbor->request(this,instance_me);
  neighbor->requests[irequest]->half = 0;
  neighbor->requests[irequest]->full = 1;
}

/* ----------------------------------------------------------------------
   init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */

double PairVMM::init_one(int i, int j)
{
  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");

  return cutmax;
}

/* ---------------------------------------------------------------------- */

void PairVMM::read_file(char *file)
{
  int params_per_line = 35;
  char **words = new char*[params_per_line+1];

  memory->sfree(params);
  params = NULL;
  nparams = maxparam = 0;

  // open file on proc 0

  FILE *fp;
  if (comm->me == 0) {
    fp = force->open_potential(file);
    if (fp == NULL) {
      char str[128];
      sprintf(str,"Cannot open Valence Multipole Model potential file %s",file);
      error->one(FLERR,str);
    }
  }

  // read each set of params from potential file
  // one set of params can span multiple lines
  // store params if all 3 element tags are in element list

  int n,nwords,ielement,jelement,kelement;
  char line[MAXLINE],*ptr;
  int eof = 0;

  while (1) {
    if (comm->me == 0) {
      ptr = fgets(line,MAXLINE,fp);
      if (ptr == NULL) {
        eof = 1;
        fclose(fp);
      } else n = strlen(line) + 1;
    }
    MPI_Bcast(&eof,1,MPI_INT,0,world);
    if (eof) break;
    MPI_Bcast(&n,1,MPI_INT,0,world);
    MPI_Bcast(line,n,MPI_CHAR,0,world);

    // strip comment, skip line if blank

    if ((ptr = strchr(line,'#'))) *ptr = '\0';
    nwords = atom->count_words(line);
    if (nwords == 0) continue;

    // concatenate additional lines until have params_per_line words

    while (nwords < params_per_line) {
      n = strlen(line);
      if (comm->me == 0) {
        ptr = fgets(&line[n],MAXLINE-n,fp);
        if (ptr == NULL) {
          eof = 1;
          fclose(fp);
        } else n = strlen(line) + 1;
      }
      MPI_Bcast(&eof,1,MPI_INT,0,world);
      if (eof) break;
      MPI_Bcast(&n,1,MPI_INT,0,world);
      MPI_Bcast(line,n,MPI_CHAR,0,world);
      if ((ptr = strchr(line,'#'))) *ptr = '\0';
      nwords = atom->count_words(line);
    }

    if (nwords != params_per_line)
      error->all(FLERR,"Incorrect format in Valence Multipole Model potential file");

    // words = ptrs to all words in line

    nwords = 0;
    words[nwords++] = strtok(line," \t\n\r\f");
    while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;

    // ielement,jelement = 1st args
    // if both args are in element list, then parse this line
    // else skip to next entry in file

    for (ielement = 0; ielement < nelements; ielement++)
      if (strcmp(words[0],elements[ielement]) == 0) break;
    if (ielement == nelements) continue;
    for (jelement = 0; jelement < nelements; jelement++)
      if (strcmp(words[1],elements[jelement]) == 0) break;
    if (jelement == nelements) continue;


    // load up parameter settings and error check their values

    if (nparams == maxparam) {
      maxparam += DELTA;
      params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
                                          "pair:params");
    }
    
    // Element_Symbol1	Element_Symbol2	We	W1	R01	B1	R02	B2	Rcut	cC1	cC2	force	De1.2	De1	De2	De3	De4	DeN	N
    params[nparams].ielement = ielement;
    params[nparams].jelement = jelement;
    strcpy(params[nparams].isymbol, words[0]);
    strcpy(params[nparams].jsymbol, words[1]);
    params[nparams].ideal_valence = atof(words[2]);
    params[nparams].k_vv = atof(words[3]);
    params[nparams].ideal_vv = atof(words[4]);
    params[nparams].w1 = atof(words[5]);
    params[nparams].r01 = atof(words[6]);
    params[nparams].b1 = atof(words[7]);
    params[nparams].r02 = atof(words[8]);
    params[nparams].b2 = atof(words[9]);
    params[nparams].cut = atof(words[10]);
    params[nparams].c1 = atof(words[11]);   // move this to setup
    params[nparams].c2 = atof(words[12]);   //move this to setup
    params[nparams].force_single = atof(words[13]);
    params[nparams].force_double = atof(words[14]);
    params[nparams].force_intercept = atof(words[15]);
    params[nparams].force_switch = atof(words[16]);  
    params[nparams].de_h1 = atof(words[17]);
    params[nparams].de_k1 = atof(words[18]);
    params[nparams].de_h2 = atof(words[19]);
    params[nparams].de_k2 = atof(words[20]);
    params[nparams].de_s02 = atof(words[21]);
    params[nparams].de_h3 = atof(words[22]);
    params[nparams].de_k3 = atof(words[23]);
    params[nparams].de_s03 = atof(words[24]);
    params[nparams].rmin_a1 = atof(words[25]);
    params[nparams].rmin_b1 = atof(words[26]);
    params[nparams].rmin_c1 = atof(words[27]);
    params[nparams].rmin_a2 = atof(words[28]);
    params[nparams].rmin_b2 = atof(words[29]);
    params[nparams].rmin_c2 = atof(words[30]);
    params[nparams].rmin_a3 = atof(words[31]);
    params[nparams].rmin_b3 = atof(words[32]);
    params[nparams].rmin_c3 = atof(words[33]);
    params[nparams].rmin_c0 = atof(words[34]);
    

    // its reall important the b values not be zero even if there is no BV term for that pair. 
    if (params[nparams].ideal_valence < 0.0 || params[nparams].k_vv < 0.0 ||
        params[nparams].ideal_vv < 0.0 || 
        params[nparams].w1 < 0.0 || params[nparams].w1 > 1.0 ||
        params[nparams].r01 < 0.0 || params[nparams].b1 < 0.0 ||
        params[nparams].r02 < 0.0 || params[nparams].b2 < 0.0 ||
        params[nparams].cut < 0.0 || params[nparams].c1 < 0.0 ||
        params[nparams].c2 < 0.0 || params[nparams].force_single < 0.0 ||
        params[nparams].de_h1 < 0.0 || params[nparams].de_k1 < 0.0 ||
        params[nparams].de_h2 < 0.0 || params[nparams].de_k2 < 0.0 ||
        params[nparams].de_k3 < 0.0 || params[nparams].de_s02 < 0.0 ||
        params[nparams].de_s03 < 0.0 || params[nparams].rmin_a1 < 0.0 ||
        params[nparams].rmin_b1 < 0.0 ||params[nparams].rmin_c1 < 0.0 ||
        params[nparams].rmin_a2 < 0.0 ||params[nparams].rmin_b2 < 0.0 ||
        params[nparams].rmin_c2 < 0.0 ||params[nparams].rmin_a3 < 0.0 ||
        params[nparams].rmin_b3 < 0.0 ||params[nparams].rmin_c3 < 0.0 ||
        params[nparams].rmin_c0 < 0.0 )

      error->all(FLERR,"Illegal Valence Multipole Model parameter");
  
    nparams++;
  }

  delete [] words;
}

/* ---------------------------------------------------------------------- */

void PairVMM::setup()
{
  int i,j,k,m,n;
  double rtmp;

  // check for duplicate entries

  memory->destroy(elem2param);
  memory->create(elem2param,nelements,nelements,"pair:elem2param");

  for (i = 0; i < nelements; i++)
    for (j = 0; j < nelements; j++) {
        n = -1;
        for (m = 0; m < nparams; m++) {
          if (i == params[m].ielement && j == params[m].jelement) {
            if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
            n = m;
          }
        }
        if (n < 0) error->all(FLERR,"Potential file is missing an entry");
        elem2param[i][j] = n;
      }

  // Set up element flags so we keep track of atomic numbers. Will be necessary to correctly compute VVS.

  for (m = 0; m < nparams; m++) {
    params[m].ielflag=0;
    params[m].jelflag=0;
    
    if (strcmp(params[m].isymbol,"O")==0) {
      params[m].ielflag=8;
    }
    if (strcmp(params[m].isymbol,"Ox")==0) {
      params[m].ielflag=108;
    }
    if (strcmp(params[m].isymbol,"H")==0) {
      params[m].ielflag=1;
    }    
    if (strcmp(params[m].isymbol,"Al")==0) {
      params[m].ielflag=13;
    } 
    if (strcmp(params[m].isymbol,"Si")==0) {
      params[m].ielflag=14;
    }    
    
    if (strcmp(params[m].jsymbol,"O")==0) {
      params[m].jelflag=8;
    }
    if (strcmp(params[m].jsymbol,"Ox")==0) {
      params[m].jelflag=108;
    }
    if (strcmp(params[m].jsymbol,"H")==0) {
      params[m].jelflag=1;
    }    
    if (strcmp(params[m].jsymbol,"Al")==0) {
      params[m].jelflag=13;
    } 
    if (strcmp(params[m].jsymbol,"Si")==0) {
      params[m].jelflag=14;
    } 


  }  

  // set cutmax to max of all params. Lammps uses this to set up neighbor lists.

  cutmax = 0.0;
  for (m = 0; m < nparams; m++) {
    if (params[m].force_switch >= 0) {
      rtmp = 2*params[m].cut;
      if (rtmp > cutmax) cutmax = rtmp;
    }
    else {
      rtmp = params[m].cut;
      if (rtmp > cutmax) cutmax = rtmp;    
    }
  }
}

/* ---------------------------------------------------------------------- */

double PairVMM::bondvalence(Param *param, double r)
{
  double sij;  
  sij=0.000000000000000000000;

  
//printf ("rcut = %f \n",param->cut);
//printf ("w1 = %f \n",param->w1);
//printf ("r01 = %f \n",param->r01);
//printf ("b1 = %f \n",param->b1);
//printf ("c1 = %f \n",param->c1);
//printf ("r02 = %f \n",param->r02);
//printf ("b2 = %f \n",param->b2);
//printf ("c2 = %f \n",param->c2);
  
  // See manual page for explanation.  
  if (r < param->cut) {
    if (param->b1 >0.0) {
      sij = param->w1*exp((param->r01-r)/param->b1)-param->c1;
    }
     if (param->b2 >0.0) {
       sij += (1-param->w1)*exp((param->r02-r)/param->b2)-param->c2;
     }
// printf ("sij1 = %f \n",sij1);
//  printf ("sij2 = %f \n",sij2);

  
  if (sij < 0) sij = 0.000000000000000000;
  
  }
  return sij;
}

/* ---------------------------------------------------------------------- */

//double PairVMM::bvgrad(double r01, double r02, double b1, double b2, double w1, double r)
//{
//  double dsij;  
//  dsij=0;

  
//  printf ("w1 = %f \n",w1);
//  printf ("r01 = %f \n",r01);
//  printf ("b1 = %f \n",b1);
//  printf ("r02 = %f \n",r02);
//  printf ("b2 = %f \n",b2);
  
  // See manual page for explanation.  
  
//    if (b1 >0.0) {
//      dsij = -(w1/b1)*exp((r01-r)/b1);
//    }
//     if (b2 >0.0) {
//       dsij += -((1-w1)/b2)*exp((r02-r)/b2);
//     }
//  printf ("dsij = %f \n",dsij);

  
//  if (dsij > 0) dsij = 0.000000000000000000;
  
  
//  return dsij;
//}
/* ---------------------------------------------------------------------- */

void PairVMM::vectorvalence(Param *param, int jnum, double *ssort, double *sjk, int jj, int *jlist, int ncoordi, 
                             double xtmp, double ytmp, double ztmp, int *jthelement, double stotal, 
                             double &vvsumj, double &vvsideal, double &kvvs) 
{
  int kk, j, k;
  double r,delx,dely,delz,vvx,vvy,vvz;
  double **x = atom->x;

  vvsumj=0.0;
  vvsideal = 0.0;
  kvvs = 0.0;
  vvx=0;
  vvy=0;
  vvz=0;

/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */    
// calculate the actual vector valence
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */

  if (stotal > 0) {
    vvsideal = param->ideal_vv;
    kvvs = param->k_vv;
    for (kk = 0; kk < jnum; kk++) {
      j = jlist[jj];
      j &= NEIGHMASK;
      k = jlist[kk];
      k &= NEIGHMASK;
      if (kk==jj) {
        delx = xtmp - x[j][0];
        dely = ytmp - x[j][1];
        delz = ztmp - x[j][2];
      }
      else {
        delx = x[k][0] - x[j][0];
        dely = x[k][1] - x[j][1];
        delz = x[k][2] - x[j][2];
      }
      r = sqrt(delx*delx + dely*dely + delz*delz);

      // k = jlist[kk];
      // k &= NEIGHMASK;
      if (sjk[kk] > 0) {

        // Here we take the bond valence to calculate the vector valence.
        // jthelement[kk]=param->jelflag;      
        vvx += delx*sjk[kk]/r;
        vvy += dely*sjk[kk]/r;
        vvz += delz*sjk[kk]/r;

        // normalize the totals to the bond valence sum.    
        vvsumj = sqrt(vvx*vvx+vvy*vvy+vvz*vvz)/stotal;  

      }
    }
  
      
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */    
// calculate the idea vector valence and force constant
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */

    if (param->ielflag == 8) {
      // This is all very primitive. Unfortuantely at the moment we have to calcualte a integer
      // coordination number. in future we will need to calcualte a real number.    
      ovvsideal(jnum, ncoordi, kvvs, vvsideal, stotal, ssort, jthelement);        
    }
    // else  {
      // This is more straightforward if H is bound primarily to oxygen it has a vvs if it isn't it doesn't.
      // if (param->ielflag == 1) {
      //   if (maxs(jnum,ssort,jthelement) != 8 )   {
      //     vvsideal = 0;
      //     kvvs = 0;
      //   } 
      //   else { 
      //     hvvsideal(jnum, ncoordi, kvvs, vvsideal, stotal, ssort, jthelement);  
      //   }
      // }
    // }
  }
}

/* ---------------------------------------------------------------------- */

void PairVMM::nonbondedenergy(Param *param, double r, int eflag, int is13bonded, double &eng,  
                      double delx, double dely, double delz, double &forcex, double &forcey, double &forcez,int idebug)
{

  double f2di, f2dj, smini, sminj, debondi, debondj, fbondi,fbondj, rmini, rminj; 
  double deij_drij, deji_drij, drij_dxi, drij_dyi, drij_dzi, devdw, debondi0, debondj0;
  int i,j,jj,nn;
  double bondenergy, rshift, xshift, bondenergy_shifti, bondenergy_shiftj, deij_drij_shift, deji_drij_shift, ai, bi, aj, bj;
  double qqrd2e = force->qqrd2e;
  double *special_coul = force->special_coul;
  double forcecoul,fpaircoul,ecoul,forcecouli,fpaircouli,ecouli,forcecoulj,fpaircoulj,ecoulj,factor_coul;


  eng=0;
  f2di=0;
  f2dj=0;
  debondi=0;
  debondj=0;
  fbondi=0;
  fbondj=0;
  rmini=0;
  rminj=0;
  deij_drij=0;
  deji_drij=0;
  devdw=0;
  bondenergy=0;
  ecoul=0;
  fpaircoul=0;

  //printf ("Inside two body. \n");

  delx=2*delx/r;
  dely=2*dely/r;
  delz=2*delz/r;
  
//  printf("sbond = %f \n",sbond);
  if (r <= param->cut) 
  {

    nonbondedparameters(rmini, debondi, fbondi, param->ielflag, param->jelflag,is13bonded);
    //nonbondedparameters(rminj, debondj, fbondj, param->jelflag, param->ielflag,is13bonded); 
    
    if(idebug ==1) {
      printf("inside non-bonded energy \n");
      printf("rmini = %f \n",rmini);
      printf("debondi = %f \n",debondi);
      printf("fbondi = %f \n",fbondi);
      printf("param->ielflag = %i \n",param->ielflag);
      printf("param->jelflag = %i \n",param->jelflag);
      printf("r = %f \n",r);
    }
    if (r > rmini && is13bonded ==1) debondi=0;
    //if (rminj > r && is13bonded ==1) debondj=0;
    
    if (debondi >0) f2di = sqrt(fbondi/(2*debondi));
    //if (debondj >0) f2dj = sqrt(fbondj/(2*debondj));

    if (is13bonded ==1) {
    // need to add charges to remove them. 
// Only owrks for a dielectric constant of 1 (vacuum)
    // factor_coul = special_coul[sbmask(j)];
    // forcecoul = qqrd2e *atom->q[i]*atom->q[j]/r;
    // fpaircoul = factor_coul*forcecoul /(r*r);
    // ecoul = factor_coul * atom->q[i]*atom->q[j]/r;
    }

    bondenergy=debondi*(exp(2*f2di*(rmini-r))-2*exp(f2di*(rmini-r))) + is13bonded*debondi;
                 //  + 0.5*debondj*(exp(2*f2dj*(rminj-r))-2*exp(f2dj*(rminj-r))) - is13bonded*debondj;  

    deij_drij=f2di*debondi*(exp(2*f2di*(rmini-r))-exp(f2di*(rmini-r)));
    //deji_drij=f2dj*debondj*(exp(2*f2dj*(rminj-r))-exp(f2dj*(rminj-r)));

    if(idebug ==1) {
      printf("bondenergy = %f \n",bondenergy);
      printf("debondi = %f \n",debondi);
      printf("is13bonded = %i \n",is13bonded);
      printf("deij_drij = %f \n",deij_drij);
    }

    // this creates a 10% smoothing zone in the outer
    // reagons of the bond. 
    rshift=0.9*param->cut;

    if (r > rshift && is13bonded==0) {
          if (is13bonded ==1) {   // This makes no sense to get to this line is13bonded=0.
    // need to add charges to remove them. 
// Only owrks for a dielectric constant of 1 (vacuum)
    factor_coul = special_coul[sbmask(j)];
    forcecoul = qqrd2e * atom->q[i]*atom->q[j]/r;
    fpaircoul = factor_coul*forcecoul /(r*r);
    ecoul = qqrd2e * factor_coul * atom->q[i]*atom->q[j]/r;
    }

      bondenergy_shifti=debondi*(exp(2*f2di*(rmini-rshift))-2*exp(f2di*(rmini-rshift)))  +is13bonded*ecoul;
      //bondenergy_shiftj=0.5*debondj*(exp(2*f2dj*(rminj-rshift))-2*exp(f2dj*(rminj-rshift))) ; 
      deij_drij_shift=f2di*debondi*(exp(2*f2di*(rmini-rshift))-exp(f2di*(rmini-rshift))) +is13bonded*fpaircoul;
      //deji_drij_shift=f2dj*debondj*(exp(2*f2dj*(rminj-rshift))-exp(f2dj*(rminj-rshift)));
                            
      xshift=param->cut-rshift;
                
      ai=(3*bondenergy_shifti-deij_drij_shift*xshift)/(xshift*xshift);
      bi=(deij_drij_shift*xshift-2*bondenergy_shifti)/(xshift*xshift*xshift);

      //aj=(3*bondenergy_shiftj-deji_drij_shift*xshift)/(xshift*xshift);
      //bj=(deji_drij_shift*xshift-2*bondenergy_shiftj)/(xshift*xshift*xshift);
                            
      bondenergy_shifti=(ai*((param->cut-r)*(param->cut-r))+bi*((param->cut-r)*(param->cut-r)*(param->cut-r)));
      deij_drij=(2*ai*(param->cut-r)+3*bi*((param->cut-r)*(param->cut-r)));
      
      //bondenergy_shiftj=(aj*((param->cut-r)*(param->cut-r))+bj*((param->cut-r)*(param->cut-r)*(param->cut-r)));
      //deji_drij=(2*aj*(param->cut-r)+3*bj*((param->cut-r)*(param->cut-r))); 
      
      bondenergy=bondenergy_shifti; //+bondenergy_shiftj;
    }     
    

    drij_dxi=delx;
    drij_dyi=dely;
    drij_dzi=delz;
      
    forcex=(deij_drij)*drij_dxi;  //+deji_drij
    forcey=(deij_drij)*drij_dyi;  //+deji_drij
    forcez=(deij_drij)*drij_dzi;  //+deji_drij
                       
          
    if (eflag) eng += bondenergy;
  }
}


/* ---------------------------------------------------------------------- */

void PairVMM::valencemonopole(Param *param, double ideal_valence_j, double r, int jnum,
                      int eflag, int is13bonded, double &eng,  double stotali, double stotalj, double sbond, 
                      double *ssorti, double *ssortj, double *sideali, double *sidealj, int ncoordi, int ncoordj,
                      double *delxj, double *delyj, double *delzj,
                      double delx, double dely, double delz, double &forcex, double &forcey, double &forcez, int i, int j, int idebug) 
{

  double f2di, f2dj, smini, sminj, debondi, debondj, fbondi,fbondj, rmini, rminj; 
  double deij_drij, deji_drij, drij_dxi, drij_dyi, drij_dzi, isnonbond;
  int jj,nn;
  double bondenergy, rshift, xshift, bondenergy_shifti, bondenergy_shiftj, deij_drij_shift, deji_drij_shift, ai, bi, aj, bj;
  double qqrd2e = force->qqrd2e;
  double *special_coul = force->special_coul;
  double forcecoul,fpaircoul,ecoul,forcecouli,fpaircouli,ecouli,forcecoulj,fpaircoulj,ecoulj,factor_coul;

  smini=0;
  sminj=0;
  f2di=0;
  f2dj=0;
  debondi=0;
  debondj=0;
  fbondi=0;
  fbondj=0;
  rmini=0;
  rminj=0;
  deij_drij=0;
  deji_drij=0;
  bondenergy=0;
  isnonbond=0;
  ecouli=0;
  ecoulj=0;
  fpaircouli=0;
  fpaircoulj=0;
  ecoul=0;
  fpaircoul=0;
  
  if (idebug==1) {
    printf ("Inside two body. \n");
    printf("sbond = %f \n",sbond);

  }

  delx=2*delx/r;
  dely=2*dely/r;
  delz=2*delz/r;
  

  if (r <= param->cut) 
  {
    for (nn = 0; nn < jnum; nn++) {

      if (sbond == ssorti[nn]) {
        smini = param->ideal_valence*sideali[nn]; 
        if (idebug==1) {
          printf("sssorti[nn] = %f \n",ssorti[nn]);
          printf("sidealitwobody[nn] = %f \n",sideali[nn]); 
        }
      }
     
      if (sbond == ssortj[nn]) {
        sminj = ideal_valence_j*sidealj[nn]; 
        if (idebug==1) {
          printf("sssortj[nn] = %f \n",ssortj[nn]);
          printf("sidealjtwobody[nn] = %f \n",sidealj[nn]); 
        }
      }
    }

    if (smini > 0.0)
    {

      if (param->w1 == 1)  
      {
        rmini = param->r01-param->b1*log(smini+param->c1);
      }
      else 
      {
        rmini =param->rmin_a1*exp((param->rmin_b1-smini)/param->rmin_c1) + param->rmin_a2*exp((param->rmin_b2-smini)/param->rmin_c2) + 
                param->rmin_a3*exp((param->rmin_b3-smini)/param->rmin_c3) + param->rmin_c0;
      } 

      // debondi0 = param->de_h1/2+param->de_h2/(1+exp(-param->de_k2*(0-param->de_s02))) +
      //    param->de_h3/(1+exp(-param->de_k3*(0-param->de_s03))) ;    
      
      debondi = param->de_h1/(1+exp(-param->de_k1*(smini)))+param->de_h2/(1+exp(-param->de_k2*(smini-param->de_s02))) +
         param->de_h3/(1+exp(-param->de_k3*(smini-param->de_s03)));    
      
      if ((smini >= param->force_switch) && (param->force_switch > 0)) {
        fbondi = (param->force_double*(smini-param->force_switch)+param->force_intercept);
      }
      else {
        fbondi = param->force_single*smini;
      }
    }
    else {
      nonbondedparameters(rmini, debondi, fbondi, param->ielflag, param->jelflag,is13bonded);
      isnonbond=1;
      if (r > rmini) {debondi=0; fbondi=0;}
    }

    
    
    if (sminj>0.0) {
      if (param->w1 == 1)   {
        rminj = param->r01-param->b1*log(sminj+param->c1);
      }
      else  {
        rminj =param->rmin_a1*exp((param->rmin_b1-sminj)/param->rmin_c1) + param->rmin_a2*exp((param->rmin_b2-sminj)/param->rmin_c2) + 
                param->rmin_a3*exp((param->rmin_b3-sminj)/param->rmin_c3) + param->rmin_c0;       
      } 

      // debondj0 = param->de_h1/2+param->de_h2/(1+exp(-param->de_k2*(0-param->de_s02))) +
      //    param->de_h3/(1+exp(-param->de_k3*(0-param->de_s03))) ;    
      
      debondj = param->de_h1/(1+exp(-param->de_k1*(sminj)))+param->de_h2/(1+exp(-param->de_k2*(sminj-param->de_s02))) +
         param->de_h3/(1+exp(-param->de_k3*(sminj-param->de_s03))) ;

      
      if ((sminj >= param->force_switch) && (param->force_switch > 0)) {
        fbondj = (param->force_double*(sminj-param->force_switch)+param->force_intercept);
      }
      else {
        fbondj = param->force_single*sminj;
      }    
    }
    else {
      nonbondedparameters(rminj, debondj, fbondj, param->jelflag, param->ielflag,is13bonded); 
      isnonbond=1;
      if (r > rminj) {debondj=0; fbondj=0;}
    }

// need to add charges to remove them. 
// Only owrks for a dielectric constant of 1 (vacuum)
    factor_coul = special_coul[sbmask(j)];
    forcecoul = qqrd2e * atom->q[i]*atom->q[j]/r;
    fpaircoul = factor_coul*forcecoul /(r*r);
    ecoul = factor_coul * qqrd2e * atom->q[i]*atom->q[j]/r;


    if (debondi >0) f2di = sqrt(fbondi/(2*debondi));
    if (debondj >0) f2dj = sqrt(fbondj/(2*debondj));


    bondenergy=0.5*debondi*(exp(2*f2di*(rmini-r))-2*exp(f2di*(rmini-r))) 
                   + 0.5*debondj*(exp(2*f2dj*(rminj-r))-2*exp(f2dj*(rminj-r))) -(1-isnonbond)*ecoul;   

    deij_drij=0.5*f2di*debondi*(exp(2*f2di*(rmini-r))-exp(f2di*(rmini-r))) -(1-isnonbond)*fpaircoul;
    deji_drij=0.5*f2dj*debondj*(exp(2*f2dj*(rminj-r))-exp(f2dj*(rminj-r))) -(1-isnonbond)*fpaircoul;


    // this creates a 10% smoothing zone in the outer
    // reagons of the bond. 
    rshift=0.9*param->cut;
    if (idebug==1) {
      printf("bondenergy = %f \n",bondenergy);
      printf("debondi = %f \n",debondi);
      printf("debondj = %f \n",debondj);
      printf("fbondi = %f \n",fbondi);
      printf("fbondj = %f \n",fbondj);
      printf("r = %f \n",r);
      printf("rmini = %f \n",rmini);
      printf("rminj = %f \n",rminj);
      printf("smini = %f \n",smini);
      printf("sminj = %f \n",sminj);
      printf("rshift = %f \n",rshift);
      printf("forcecoul = %f \n",forcecoul);
      printf("factorcoul = %f \n",factor_coul);
      printf("q(i) = %f \n",atom->q[i]);
      printf("q(j) = %f \n",atom->q[j]);
      printf("qqrd2e = %f \n",qqrd2e);
      printf("ecoul = %f \n",ecoul);
      printf("fpaircoul = %f \n",fpaircoul);
    }

// need to add charge here too

    if (r > rshift )  /// && smini >0 && sminj > 0) This allows non-bonded zero bonds to be smoothed too. Its not ideal but untill we can generate a long range term its the best for now. 
    {

      // first charge correction
      forcecouli = qqrd2e *atom->q[i]*atom->q[j]/(rmini-rshift);
      fpaircouli = factor_coul*forcecouli /((rmini-rshift)*(rmini-rshift));
      ecouli = factor_coul * atom->q[i]*atom->q[j]/(rmini-rshift);
      forcecoulj = qqrd2e * atom->q[i]*atom->q[j]/(rminj-rshift);
      fpaircoulj = factor_coul*forcecoulj /((rminj-rshift)*(rminj-rshift));
      ecoulj = factor_coul * qqrd2e * atom->q[i]*atom->q[j]/(rminj-rshift);
    
      // then inital boundary conditions far boundary zero
      bondenergy_shifti=0.5*debondi*(exp(2*f2di*(rmini-rshift))-2*exp(f2di*(rmini-rshift))) -0.5*(1-isnonbond)*ecouli;
      bondenergy_shiftj=0.5*debondj*(exp(2*f2dj*(rminj-rshift))-2*exp(f2dj*(rminj-rshift))) -0.5*(1-isnonbond)*ecoulj; 
      deij_drij_shift=0.5*f2di*debondi*(exp(2*f2di*(rmini-rshift))-exp(f2di*(rmini-rshift)))-(1-isnonbond)*fpaircouli;
      deji_drij_shift=0.5*f2dj*debondj*(exp(2*f2dj*(rminj-rshift))-exp(f2dj*(rminj-rshift)))-(1-isnonbond)*fpaircoulj;
                            
      xshift=param->cut-rshift;
                
      ai=(3*bondenergy_shifti-deij_drij_shift*xshift)/(xshift*xshift);
      bi=(deij_drij_shift*xshift-2*bondenergy_shifti)/(xshift*xshift*xshift);

      aj=(3*bondenergy_shiftj-deji_drij_shift*xshift)/(xshift*xshift);
      bj=(deji_drij_shift*xshift-2*bondenergy_shiftj)/(xshift*xshift*xshift);
                            
      bondenergy_shifti=(ai*((param->cut-r)*(param->cut-r))+bi*((param->cut-r)*(param->cut-r)*(param->cut-r)));
      deij_drij=(2*ai*(param->cut-r)+3*bi*((param->cut-r)*(param->cut-r)));
      
      bondenergy_shiftj=(aj*((param->cut-r)*(param->cut-r))+bj*((param->cut-r)*(param->cut-r)*(param->cut-r)));
      deji_drij=(2*aj*(param->cut-r)+3*bj*((param->cut-r)*(param->cut-r))); 
      
      bondenergy=bondenergy_shifti+bondenergy_shiftj;
    }     

    drij_dxi=delx;
    drij_dyi=dely;
    drij_dzi=delz;
      
    forcex=(deij_drij+deji_drij)*drij_dxi;  
    forcey=(deij_drij+deji_drij)*drij_dyi;
    forcez=(deij_drij+deji_drij)*drij_dzi; 
                       
          
    if (eflag) eng += bondenergy;
    
  }
}

/* ---------------------------------------------------------------------- */

void PairVMM::valencedipole(Param *param, double r,
                      int eflag, double &eng,  double stotal, double sbond,
                      double vvstotal, double vvsideal, double vvx, double vvy, double vvz, double kvvs,
                      double delx, double dely, double delz, double &forcex, double &forcey, double &forcez, int idebug)
{
//  double rinvsq,rp,rq,rainv,rainvsq,expsrainv;
  double fmagnitude, vnorm, vdotb,bdotb;
  double vfx, vfy, vfz;
  
// printf ("Inside two body. \n");

// printf("3B vvsideal = %f \n",vvsideal);
// printf("3B kvvs = %f \n",kvvs);
// printf("2B vvstotal = %f \n",vvstotal);

  delx=2*delx/r;
  dely=2*dely/r;
  delz=2*delz/r;
 
  if (r <= param->cut) 
    {
            
    fmagnitude =     (sbond/stotal)*0.5*kvvs*(vvstotal-vvsideal);
    vdotb = vvx*delx+vvy*dely+vvz*delz;
    bdotb = delx*delx+dely*dely+delz*delz;
      
//      printf("fmagnitude = %f \n",fmagnitude);
//      printf("vdotb = %f \n",vdotb);
//      printf("bdotb = %f \n",bdotb);
    vfx = vvx - (vdotb/bdotb)*delx;
    vfy = vvy - (vdotb/bdotb)*dely;
    vfz = vvz - (vdotb/bdotb)*delz;
    vnorm = sqrt(vfx*vfx+vfy*vfy+vfz*vfz);
//      printf("2B vnorm = %f \n",vnorm);
    if (vnorm > 0)
      {
      forcex=-fmagnitude*(vfx)/vnorm;
      forcey=-fmagnitude*(vfy)/vnorm;
      forcez=-fmagnitude*(vfz)/vnorm;
      }
    else
      {
        forcex=0;
        forcey=0;
        forcez=0;
      }
    }  
      
//      printf("VV force x = %f \n",forcex);   
//      printf("VV force y = %f \n",forcey);   
//      printf("VV force z = %f \n",forcez); 
          
    if (eflag) eng += kvvs*(vvstotal-vvsideal)*(vvstotal-vvsideal);
    if (idebug ==1) {
      printf("Valence Vector energy = %f \n",kvvs*(vvstotal-vvsideal)*(vvstotal-vvsideal));
    }

    
}

/* ---------------------------------------------------------------------- */

int PairVMM::ncoord(int jnum, double *sj, double stotal, int *jthelement)
{
  int jmax,ii,jj, ncord;
  double stemp;
 
  jmax=jnum;
  ncord=0;
  stemp=0;
  
  // This code not only computes the coordination number as an integer it also computes the strongest bonds, the elements atomic numbers
  // corresponding to that type. 
  if (jnum >0)
  {
  
    for (ii = 0; ii < jnum; ii++) {
      stemp+=sj[ii]*sj[ii]/(stotal*stotal); 
      // printf("Sj[ii] (before) = %f \n",sj[ii]);
      // printf("stotal = %f \n",stotal);
      // printf("Stemp (accumulating) = %f \n",stemp);
      
    }
    if (stemp > 0) {
      stemp=1/(stemp);
    // printf("Stemp = %f \n",stemp);
      ncord=int(ceil(stemp-0.02));
    // printf("ncord = %i \n",ncord);
    }
       
    for (ii = 0; ii < jnum-1; ii++) {
      for (jj = ii+1; jj < jnum; jj++) {

        if (sj[ii] < sj[jj]) {
          std::swap(sj[ii],sj[jj]);
          std::swap(jthelement[ii],jthelement[jj]);
        
        }
          //  printf("Sj[ii] (after) = %f \n",sj[ii]);

      }
    }
  }
  return ncord; 
} 

/* ---------------------------------------------------------------------- */

int PairVMM::maxs(int jnum, double *sj, int *jthelement) 
{
  // This computes the atomic number of the bond corresponding to the strongest bond. Used for hydrogen.
  double smax;
  int jmax,jj;
  
  smax=0;
  jmax=0;
  
  for (jj = 0; jj < jnum; jj++) {
    if (sj[jj] > smax) {
      smax = sj[jj];
      jmax = jthelement[jj]; 

    }
  }
  return jmax;     
} 



/* ---------------------------------------------------------------------- */

void PairVMM::ovvsideal(int jnum, int ncord, double &kvvs, double &vvsideal, double stot, double *ssort, int *jthelement)
                       
{
  // This code computes the ideal vvs for oxygen based on what it is bound to.
  double angle1, angle2, angle3, PI;
  int jsum;
  
  PI=3.14159265359;
  jsum=0;
  
//  printf("ncord = %i \n",ncord);
  if (ncord ==1)
  {
    kvvs=0;
    vvsideal=0; 
  }
  else if (ncord ==2)
  {
      jsum = (jthelement[0] + jthelement[1]) % 100;
 //     printf("jsum = %i \n",jsum);
      if (jsum == 2)  angle1=104.5;
      if (jsum == 9) angle1=101.9;
      if (jsum == 14) angle1=130.8;
      if (jsum == 15) angle1=118.7;
      if (jsum == 16) angle1=116.78;
      if (jsum == 21) angle1=109;  // fix
      if (jsum == 22) angle1=140; // fix
      if (jsum == 26) angle1=180.0;
      if (jsum == 27) angle1=180.0;
      if (jsum == 28) angle1=145.5;
      

//      kvvs=param->k_vv;
//      printf("O2 ssort[0] = %f \n",ssort[0]);
//      printf("O2 ssort[1] = %f \n",ssort[1]);

      vvsideal = sqrt(pow(ssort[0],2)+pow(ssort[1],2) 
           +2*ssort[0]*ssort[1]*cos(angle1*PI/180.0))/stot;   
//      printf("O2 vvsideal = %f \n",vvsideal);
//      printf("O2 kvvs = %f \n",kvvs); 
//     if (ssort[1] > 0.0001) {
//       kvvs=kvvs/(ssort[0]*sqrt(ssort[1]));
//     }
  }
  else if (ncord ==3)
  {
     jsum =(jthelement[0] + jthelement[1] + jthelement[2]) % 100 ;
//     printf("jsum = %i \n",jsum);
     
     angle1=115;
     angle2=115;
     angle3=115;
       
     if (jsum ==3) {
         angle1=111.3;
         angle2=111.3;
         angle3=111.3;
         
         if (ssort[2] < 0.2)  angle1=108.7163383;
         
     }

     if ((jsum >= 39) && (jsum <=42)) {
         angle1=98.6;
         angle2=130.7;
         angle3=130.7;
     }
//     printf("angle1 = , %f \n",angle1);
//     printf("angle2 = , %f \n",angle2);
//     printf("angle3 = , %f \n",angle3);
    
     vvsideal=sqrt(pow(ssort[0],2)+pow(ssort[1],2)+pow(ssort[2],2) 
           +2*ssort[0]*ssort[1]*cos(angle1*PI/180.0)+2*ssort[0]*ssort[2]*cos(angle2*PI/180.0)
           +2*ssort[1]*ssort[2]*cos(angle3*PI/180.0))/stot;
//    if (ssort[1] > 0.0001) {
//      kvvs=kvvs/(ssort[0]*sqrt(ssort[1]));
//     }
//     printf("O3 vvsideal = %f \n",vvsideal);
//     printf("O3 kvvs = %f \n",kvvs); 
     
  }  
  else if  (ncord == 4)
  {
 //   if (ssort[1] > 0.0001) {
 //     kvvs=kvvs/(ssort[0]*sqrt(ssort[1]));
 //   }
//      vvsideal=param->ideal_vv; 
//      printf("O4 vvsideal = %f \n",vvsideal);
//      printf("O4 kvvs = %f \n",kvvs); 
  }
  else if (ncord > 4)
  {
    kvvs=0;
    vvsideal=0; 
//      printf("O5 vvsideal = %f \n",vvsideal);
//      printf("O5 kvvs = %f \n",kvvs); 
  } 
    
//      printf("O vvsideal = %f \n",vvsideal);
//      printf("O kvvs = %f \n",kvvs); 
     
} 

/* ---------------------------------------------------------------------- */

void PairVMM::hvvsideal(int jnum, int ncord, double &kvvs, double &vvsideal, double stot, double *ssort, int *jthelement)
                       
{
  // This code computes the ideal vvs for oxygen based on what it is bound to.
  double angle1, angle2, angle3, PI;
  int jsum;
  
  PI=3.14159265359;
  angle1=172.0;

//  printf("ncord = %i \n",ncord);

  jsum = (jthelement[0] + jthelement[1]) % 100;

//  if (jsum == 16) 
//    {
//      angle1=173.0;
//      if (ssort[1] > 0.001) {
//        kvvs=kvvs/(ssort[0]*sqrt(ssort[1]));
//      }
//    }
//  else
//    {
      kvvs=0;
      vvsideal=0; 
//    }

//      kvvs=param->k_vv;
//      printf("O2 ssort[0] = %f \n",ssort[0]);
//      printf("O2 ssort[1] = %f \n",ssort[1]);

      vvsideal = sqrt(pow(ssort[0],2)+pow(ssort[1],2) 
           +2*ssort[0]*ssort[1]*cos(angle1*PI/180.0))/stot;   
//      printf("O2 vvsideal = %f \n",vvsideal);
//      printf("O2 kvvs = %f \n",kvvs); 
     
} 

/* ---------------------------------------------------------------------- */

void PairVMM::idealconfiguration(int ncoord,int jnum, double valence, double *ssort, double stotal, int elflag, int * jelflag, double *sideal)
{
  int jj;
  double ctemp0,ctemp1,ctemp2,ctemp3;
  double saltconfig[1000];
  
  ctemp0=0;
  ctemp1=0;
  ctemp2=0;
  ctemp3=0;
//  printf("inside ideal configuration \n");
//  printf("stotal = %f \n ",stotal);
//  printf("ncoord = %i \n ",ncoord);
  

  for (jj = 0; jj < jnum; jj++) {
    sideal[jj] = 0; 
    saltconfig[jj]=0;  
 //   printf("sidealgenerici[jj] = %f \n",sideal[jj]);    
  }

  // printf("ncoord(ic) =  %i \n",ncoord );
  if (ncoord > 0) {
    for (jj = 0; jj < ncoord; jj++) {
      sideal[jj] = 1/double(ncoord); 
 //     printf("sidealgenerici[jj] = %f \n",sideal[jj]);    
    }
  }

// Determine coordination number exceptions.
      
  for (jj = 0; jj < jnum; jj++) {
    ctemp0+=sideal[jj]*ssort[jj]/valence;
    ctemp1+=(ssort[jj]/valence)*(ssort[jj]/valence);    
  }
  ctemp0=ctemp0/sqrt(ctemp1);


  
  if ((elflag ==8  || elflag ==108 )  && (jelflag[0]==1)  && (ssort[0] > 0 ) ) {  
    saltconfig[0]=0.5;
    saltconfig[1]=0.5;
    saltconfig[2]=0.0;
    saltconfig[3]=0.0;

    // if (ctemp2 > ctemp0) {
    //   ctemp0=ctemp2;
      for (jj = 0; jj < jnum; jj++) {
        sideal[jj]=saltconfig[jj];
  //      printf("sideal[jj] = %f \n",sideal[jj]);    
      }
    }

    // for (jj = 0; jj < jnum; jj++) {
    //   ctemp2+=saltconfig[jj]*ssort[jj]/valence;
    //   ctemp3+=saltconfig[jj]*saltconfig[jj];
    //   if (jnum < 4) {
    //     for (jj = jnum; jj < 4; jj++) {
    //       ctemp3+=saltconfig[jj]*saltconfig[jj];
    //     }
    //   }
    // }
    // ctemp2=ctemp2/sqrt(ctemp3);
 //   printf("ctemp2 oo = %f \n", ctemp2);
 //   printf("ctemp0 oo = %f \n", ctemp0);
           

 // } 


// This really doesn't do anything as it is written. But neither does the matlab code so shrug.
//   if ((elflag ==8) && (jelflag ==1) && (ncoord==1 ))  
//   {
//     saltconfig[0]=0.5;
//     saltconfig[1]=0.0;
//     saltconfig[2]=0.0;
//     saltconfig[3]=0.0;
//     for (jj = 0; jj < jnum; jj++) {
//       ctemp2+=saltconfig[jj]*ssort[jj]/valence;
//       ctemp3+=saltconfig[jj]*saltconfig[jj];
//     }
//     ctemp2=ctemp2/sqrt(ctemp3);
// //   printf("ctemp2 oo = %f \n", ctemp2);
// //    printf("ctemp0 oo = %f \n", ctemp0);
           
//     if (ctemp2 > ctemp0) {
//       ctemp0=ctemp2;
//       for (jj = 0; jj < jnum; jj++) {
//         sideal[jj]=saltconfig[jj];
//       }
//     }
//   }  

  ctemp2=0;
  ctemp3=0; 
  if ((elflag ==1) && (ssort[0] > 0))
  {
    saltconfig[0]=1.00;
    saltconfig[1]=0.00;
    saltconfig[2]=0.00;
    saltconfig[3]=0.00;
//     for (jj = 0; jj < jnum; jj++) {
//       ctemp2+=saltconfig[jj]*ssort[jj]/valence;
//       ctemp3+=saltconfig[jj]*saltconfig[jj];
//     }
//     ctemp2=ctemp2/sqrt(ctemp3);
// //    printf("ctemp2 = %f \n", ctemp2);
    
           
//     if ((ctemp2 > ctemp0) || (ssort[2] < 0.2)) {
//       ctemp0=ctemp2;
      for (jj = 0; jj < jnum; jj++) {
        sideal[jj]=saltconfig[jj];
//        printf("sideal[jj] = %f \n",sideal[jj]);  
      }
    }
//  }

//   ctemp2=0;
//   ctemp3=0;
//   if ((elflag ==1) && (ncoord>2)) 
//   {
//     saltconfig[0]=0.50;
//     saltconfig[1]=0.50;
//     saltconfig[2]=0.00;
//     saltconfig[3]=0.00;
//     for (jj = 0; jj < jnum; jj++) {
//       ctemp2+=saltconfig[jj]*ssort[jj]/valence;
//       ctemp3+=saltconfig[jj]*saltconfig[jj];
//     }
//     ctemp2=ctemp2/sqrt(ctemp3);
// //    printf("ctemp2 = %f \n", ctemp2);
           
//     if (ctemp2 > ctemp0){
//       ctemp0=ctemp2;
//       for (jj = 0; jj < jnum; jj++) {
//         sideal[jj]=saltconfig[jj];
//       }
//     }
//   }
//   for (jj = 0; jj < jnum; jj++) {
//     sideal[jj]=valence*sideal[jj];
//   }
}
/* ---------------------------------------------------------------------- */

void PairVMM::nonbondedparameters(double &rmin,double &debond, double &fbond,int ielflag, int jelflag, int is13bonded) 
{
  double rmini,rminj,debondi,debondj,fbondi,fbondj,debondi0,debondj0;

  rmini=0;
  rminj=0;
  debondi=0;
  debondj=0;
  fbondi=0;
  fbondj=0;

  if (ielflag >100) {
    ielflag=ielflag-100;
  }

  if (jelflag >100) {
    jelflag=jelflag-100;
  }
  
  rmini = 0.4557*log((double)ielflag) + 2.6395;
  debondi = (0.2767*log((double)ielflag) + 0.1478)/4.184;
  fbondi = (0.8241*log((double)ielflag) + 1.4476)/4.184;
  

  rminj = 0.4557*log((double)jelflag) + 2.6395;
  debondj = (0.2767*log((double)jelflag) + 0.1478)/4.184;
  fbondj = (0.8241*log((double)jelflag) + 1.4476)/4.184;
 
  if (is13bonded==1) {
    rmini=rmini/3;
    rminj=rminj/3;
    fbondi=fbondi*3;
    fbondj=fbondj*3;
    
    // if (ielfag==1) {
    //   rmini=1.507133834;
    //   debondi=1.084553138;
    //   fbondi=24.3634582;
    //   debondi0=debondi;
    // }
    // if (jelfag==1) {
    //   rminj=1.507133834;
    //   debondj=1.084553138;
    //   fbondj=24.3634582;
    //   debondj0=debondj;
    // }
  }

  // if (ielflag==8) {
  //   rmini=3.5532;
  //   debondi=0.1554;
  //   fbondi=1.271828446;
  // }
  // if (jelflag==8) {
  //   rminj=3.5532;
  //   debondj=0.1554;
  //   fbondj=1.271828446;
  // }

  rmin=(rmini+rminj)/2;
  debond=pow(debondi*debondj,0.5);
  fbond=pow(fbondi*fbondj,0.5);
}