LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Coul/cut gradients question
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Coul/cut gradients question


From: Matthew Wander <mcfwander@...24...>
Date: Thu, 6 Jul 2017 12:51:05 -0400

Okay. Since I don't actually understand the problem it will be challenging but I have determined more of the issue. Here is what additional information I have. I am running a series of two body scans to check the gradients of my pair style. Attached is one sample input. There are actually several hundred of these and the problem exists in all of them regardless of bonds or atom types or charges or distances or anything else that I can think of. You can see clearly in the output that there is an energy contribution from the charge which there should be. However when one checks the dump file the forces on the atom are zero which is definitely not correct. I tried using hybrid overlay and two standard pair styles and that does give the correct dumping of gradients and energies. 

The only thing I can think of, and I have no idea if it is correct, is that somehow the problem is in the invocation of evtally. I use evtally_xyz_full and I know that could/cut uses the much simple could/cut. Or maybe I have a flag set correctly? All I know is that my pair style is producing the correct information for its part, but somehow is overriding the other pair style for  charge but I have no idea how that is possible.

thank you very much in advance for any ideas

matthew


On Thu, Jul 6, 2017 at 10:17 AM, Axel Kohlmeyer <akohlmey@...24...> wrote:
On Thu, Jul 6, 2017 at 10:08 AM, Matthew Wander <mcfwander@...24...> wrote:
> Hi
>
> I am using hybrid overlay to conduct a two body potential energy scan using
> pair_coul_cut and my pair style vmm. However, when I check my dump file I
> only see the forces computed by my pair style not the contribution from
> coul/cut. Coul/cut is correctly contributing to the energy but somehow the
> contributions to f[i] are not, or if they are they are not being dumped
> correctly.
>
> Any ideas?

there is not much information to go on. you are missing even the most
basic required information (e.g. LAMMPS version, platform, compiler,
example input).

have you tried running with pair style coul/cut alone?

axel.

>
> thanks
> matthew
>
>
> ------------------------------------------------------------------------------
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> _______________________________________________
> lammps-users mailing list
> lammps-users@...12...396...sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/lammps-users
>



--
Dr. Axel Kohlmeyer  akohlmey@...43...4...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.

/* ----------------------------------------------------------------------
   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);
}
/* -*- c++ -*- ----------------------------------------------------------
   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.
------------------------------------------------------------------------- */

#ifdef PAIR_CLASS

PairStyle(vmm,PairVMM)

#else

#ifndef LMP_PAIR_VMM_H
#define LMP_PAIR_VMM_H

#include "pair.h"

namespace LAMMPS_NS {

class PairVMM : public Pair {
 public:
  PairVMM(class LAMMPS *);
  virtual ~PairVMM();
  virtual void compute(int, int);
  void settings(int, char **);
  virtual void coeff(int, char **);
  virtual double init_one(int, int);
  virtual void init_style();

  struct Param {                // Element_Symbol1	Element_Symbol2	We	W1	R01	B1	R02	B2	Rcut	cC1	cC2	force	De1.2	De1	De2	De3	De4	DeN	N
    double ideal_valence,k_vv,ideal_vv;   // This row corresponds to parameter values for only the ith atom.
    double w1;
    double r01,b1,r02,b2,c1,c2;
    double force_single,force_double,force_intercept,force_switch;
    double de_h1,de_k1,de_h2,de_k2,de_s02,de_h3,de_k3,de_s03;
    double rmin_a1,rmin_b1,rmin_c1,rmin_a2,rmin_b2,rmin_c2,rmin_a3,rmin_b3,rmin_c3,rmin_c0;
    double tol;                 // remove tol later
    double cut,cutsq;           //remove custsq later
    int ielement,jelement,kelement; //remove kelement
    int ielflag,jelflag;
    char isymbol[3],jsymbol[3];
  };
  
//  struct Bond {
//    double *smin;
//    int firststep;
//  };
  
 protected:
  double cutmax;                // max cutoff for all elements
  int nelements;                // # of unique elements
  char **elements;              // names of unique elements
  int **elem2param;            // mapping from element doublets to parameters
  int *map;                     // mapping from atom types to elements
  int nparams;                  // # of stored parameter sets
  int maxparam;                 // max # of parameter sets
  Param *params;                // parameter set for an I-J-K interaction
//  Bond *bonds;                  // temp holder for valence informaiton
  
  virtual void allocate();
  void read_file(char *);
  virtual void setup();
  int ncoord(int, double *, double, int *);
  void ovvsideal(int, int, double &, double &, double, double *, int *);
  void hvvsideal(int, int, double &, double &, double, double *, int *);
  void idealconfiguration(int, int, double, double *, double, int, int *, double *);
  int maxs(int, double *, int *);  
  double bondvalence(Param *, double);
//  double bvgrad(double, double, double, double, double, double);
  void vectorvalence(Param *, int, double *, double *, int, int *,  int, double, double, double, int *, 
                      double, double &, double &, double &);
  void nonbondedenergy(Param *, double, int, int, double &, double, double, double, double &, double &, double &,int);
  void valencemonopole(Param *, double , double, int, int, int, double &,  double, double, double, 
                      double *, double *, double *, double *, int, int,
                      double *, double *, double *, 
                      double, double, double, double &, double &, double &, int, int,int);
  void valencedipole(Param *, double, int, double &, 
                  double, double, double, double, double , double , double ,double,
                  double, double, double, double &, double &, double &, int);
  void nonbondedparameters(double &,double &, double &,int,int,int); 
};

}

#endif
#endif

/* ERROR/WARNING messages:

E: Illegal ... command

Self-explanatory.  Check the input script syntax and compare to the
documentation for the command.  You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.

E: Incorrect args for pair coefficients

Self-explanatory.  Check the input script or data file.

E: Pair style Valence Multipole Model requires atom IDs

This is a requirement to use the VMM potential.

E: Pair style Valence Multipole Model requires newton pair on

See the newton command.  This is a restriction to use the VMM
potential.

E: All pair coeffs are not set

All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.

E: Cannot open Valence Multipole Model potential file %s

The specified VMM potential file cannot be opened.  Check that the path
and name are correct.

E: Incorrect format in Valence Multipole Model potential file

Incorrect number of words per line in the potential file.

E: Illegal Valence Multipole Model parameter

One or more of the coefficients defined in the potential file is
invalid.

E: Potential file has duplicate entry

The potential file has more than one entry for the same element.

E: Potential file is missing an entry

The potential file does not have a needed entry.

*/

Attachment: hh_020_vmm.data
Description: Binary data

Attachment: hh_020_vmm.inp
Description: Binary data

Attachment: hh_020_vmm.log
Description: Binary data

Attachment: hh_020_vmm.xyz
Description: Protein Databank data