void BondMyQuartic::compute(int eflag, int vflag)							//ADDED
{
  int i1,i2,n,m,type,itype,jtype;
  int ai1,ai2,ai3,am,an;				//ADDED
  double delx,dely,delz,ebond,fbond,evdwl,fpair;
  double r,rsq,dr,r2,ra,rb,sr2,sr6;

  ebond = evdwl = 0.0;
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = 0;

  // insure pair->ev_tally() will use 1-4 virial contribution

  if (vflag_global == 2)
    force->pair->vflag_either = force->pair->vflag_global = 1;

  double **cutsq = force->pair->cutsq;
  double **x = atom->x;
  double **f = atom->f;
  int **bondlist = neighbor->bondlist;
  int nbondlist = neighbor->nbondlist;
  int **anglelist = neighbor->anglelist;	//ADDED
  int nanglelist = neighbor->nanglelist;	//ADDED
  int nlocal = atom->nlocal;
  int newton_bond = force->newton_bond;

  for (n = 0; n < nbondlist; n++) {

    // skip bond if already broken

    if (bondlist[n][2] <= 0) continue;

    i1 = bondlist[n][0];
    i2 = bondlist[n][1];
    type = bondlist[n][2];

    delx = x[i1][0] - x[i2][0];
    dely = x[i1][1] - x[i2][1];
    delz = x[i1][2] - x[i2][2];
    domain->minimum_image(delx,dely,delz);

    rsq = delx*delx + dely*dely + delz*delz;

    // if bond breaks, set type to 0
    //   both in temporary bondlist and permanent bond_type
    // if this proc owns both atoms,
    //   negate bond_type twice if other atom stores it
    // if other proc owns 2nd atom, other proc will also break bond

    if (rsq > rc[type]*rc[type]) {

      bondlist[n][2] = 0;

      for (m = 0; m < atom->num_bond[i1]; m++)
	if (atom->bond_atom[i1][m] == atom->tag[i2])
	  atom->bond_type[i1][m] = 0;
      if (i2 < atom->nlocal)
	for (m = 0; m < atom->num_bond[i2]; m++)
	  if (atom->bond_atom[i2][m] == atom->tag[i1])
	    atom->bond_type[i2][m] = 0;
    //  continue;

	  /* --------------ADDED to remove angles in temporal and permanent list---------------------------------------------*/
	  for (an = 0; an < nanglelist; an++)
	  {
		  if (anglelist[an][3] <= 0) continue;

		  ai1 = anglelist[an][0];
	      	  ai2 = anglelist[an][1];
                  ai3 = anglelist[an][2];

		  if (((atom->tag[ai1] == atom->tag[i1])&&(atom->tag[ai2] == atom->tag[i2]))||((atom->tag[ai1] == atom->tag[i2])&&(atom->tag[ai2] == atom->tag[i1])))
		  {
			  anglelist[an][3] = 0;

			  if (ai1 < atom->nlocal)
			  {
				  for (am = 0; am < atom->num_angle[ai1]; am++)
				  {
					  if ((atom->angle_atom2[ai1][am] == atom->tag[ai2])||(atom->angle_atom1[ai1][am] == atom->tag[ai2])||(atom->angle_atom3[ai1][am] == atom->tag[ai2]))
					  //if ((atom->angle_atom2[ai1][am] == atom->tag[ai2])&&(atom->angle_atom3[ai1][am] == atom->tag[ai3]))
					  {
						  atom->angle_type[ai1][am] = 0;
					  }
					  //if ((atom->angle_atom2[ai1][am] == atom->tag[ai2])&&(atom->angle_atom1[ai1][am] == atom->tag[ai3]))
					  //{
						  //atom->angle_type[ai1][am] = 0;
					  //}
				  }
			  }
			  if (ai2 < atom->nlocal)
			  {
				  for (am = 0; am < atom->num_angle[ai2]; am++)
				  {
					  if ((atom->angle_atom1[ai2][am] == atom->tag[ai1])||(atom->angle_atom3[ai2][am] == atom->tag[ai1])||(atom->angle_atom2[ai2][am] == atom->tag[ai1]))
					  //if ((atom->angle_atom1[ai2][am] == atom->tag[ai1])&&(atom->angle_atom3[ai2][am] == atom->tag[ai3]))
					  {
						  atom->angle_type[ai2][am] = 0;
					  }
					  //if ((atom->angle_atom3[ai2][am] == atom->tag[ai1])&&(atom->angle_atom1[ai2][am] == atom->tag[ai3]))
					  //{
						  //atom->angle_type[ai2][am] = 0;
					  //}
				  }
			  }
			  if (ai3 < atom->nlocal)
			  {
				  for (am = 0; am < atom->num_angle[ai3]; am++)
				  {
					  if ((atom->angle_atom1[ai3][am] == atom->tag[ai1])&&(atom->angle_atom2[ai3][am] == atom->tag[ai2]))
					  {
						  atom->angle_type[ai3][am] = 0;
					  }
					  if ((atom->angle_atom3[ai3][am] == atom->tag[ai1])&&(atom->angle_atom2[ai3][am] == atom->tag[ai2]))
					  {
						  atom->angle_type[ai3][am] = 0;
					  }
					  if ((atom->angle_atom1[ai3][am] == atom->tag[ai2])&&(atom->angle_atom2[ai3][am] == atom->tag[ai1]))
					  {
						  atom->angle_type[ai3][am] = 0;
					  }
					  if ((atom->angle_atom3[ai3][am] == atom->tag[ai2])&&(atom->angle_atom2[ai3][am] == atom->tag[ai1]))
					  {
						  atom->angle_type[ai3][am] = 0;
					  }		
				  }
			  }
		  }

		  if (((atom->tag[ai2] == atom->tag[i1])&&(atom->tag[ai3] == atom->tag[i2]))||((atom->tag[ai2] == atom->tag[i2])&&(atom->tag[ai3] == atom->tag[i1])))
		  {
			  anglelist[an][3]=0;

			  if (ai1 < atom->nlocal)
			  {
				  for (am = 0; am < atom->num_angle[ai1]; am ++)
				  {
					  if ((atom->angle_atom2[ai1][am] == atom->tag[ai2])&&(atom->angle_atom3[ai1][am] == atom->tag[ai3]))
					  {
						  atom->angle_type[ai1][am] = 0;
					  }
					  if ((atom->angle_atom2[ai1][am] == atom->tag[ai2])&&(atom->angle_atom1[ai1][am] == atom->tag[ai3]))
					  {
						  atom->angle_type[ai1][am] = 0;
					  }
					  if ((atom->angle_atom2[ai1][am] == atom->tag[ai3])&&(atom->angle_atom3[ai1][am] == atom->tag[ai2]))
					  {
						  atom->angle_type[ai1][am] = 0;
					  }
					  if ((atom->angle_atom2[ai1][am] == atom->tag[ai3])&&(atom->angle_atom1[ai1][am] == atom->tag[ai2]))
					  {
						  atom->angle_type[ai1][am] = 0;
					  }		
				  }
			  }
			  if (ai2 < atom->nlocal)
			  {
				  for (am = 0; am < atom->num_angle[ai2]; am ++)
				  {
					  if ((atom->angle_atom3[ai2][am] == atom->tag[ai3])||(atom->angle_atom1[ai2][am] == atom->tag[ai3])||(atom->angle_atom2[ai2][am] == atom->tag[ai3]))
					  //if ((atom->angle_atom1[ai2][am] == atom->tag[ai1])&&(atom->angle_atom3[ai2][am] == atom->tag[ai3]))
					  {
						  atom->angle_type[ai2][am] = 0;
					  }
					  //if ((atom->angle_atom3[ai2][am] == atom->tag[ai1])&&(atom->angle_atom1[ai2][am] == atom->tag[ai3]))
					  //{
						  //atom->angle_type[ai2][am] = 0;
					  //}
				  }
			  }
			  if (ai3 < atom->nlocal)
			  {
				  for (am = 0; am < atom->num_angle[ai3]; am ++)
				  {
					  if ((atom->angle_atom2[ai3][am] == atom->tag[ai2])||(atom->angle_atom1[ai3][am] == atom->tag[ai2])||(atom->angle_atom3[ai3][am] == atom->tag[ai2]))
				  	  //if ((atom->angle_atom1[ai3][am] == atom->tag[ai1])&&(atom->angle_atom2[ai3][am] == atom->tag[ai2]))
					  {
						  atom->angle_type[ai3][am] = 0;
					  }
					  //if ((atom->angle_atom3[ai3][am] == atom->tag[ai1])&&(atom->angle_atom2[ai3][am] == atom->tag[ai2]))
					  //{
						  //atom->angle_type[ai3][am] = 0;
					  //}
				  }
			  }
		  }
	  }
	  /* --------------ADDED---------------------------------------------*/

	  continue;

    }

    // quartic bond
    // 1st portion is from quartic term
    // 2nd portion is from LJ term cut at 2^(1/6) with eps = sigma = 1.0

    r = sqrt(rsq);
    dr = r - rc[type];
    r2 = dr*dr;
    ra = dr - b1[type];
    rb = dr - b2[type];
    fbond = -k[type]/r * (r2*(ra+rb) + 2.0*dr*ra*rb);

    //if (rsq < TWO_1_3) {
	if (rsq < TWO_1_3*sigma[type]*sigma[type]){
      //sr2 = 1.0/rsq;
	  sr2 = sigma[type]*sigma[type]/rsq;
      sr6 = sr2*sr2*sr2;
      //fbond += 48.0*sr6*(sr6-0.5)/rsq;
	  fbond += 48.0*epsilon[type]*sr6*(sr6-0.5)/rsq;
    }

    if (eflag) {
      ebond = k[type]*r2*ra*rb + u0[type];
      //if (rsq < TWO_1_3) ebond += 4.0*sr6*(sr6-1.0) + 1.0;
	  if (rsq < TWO_1_3*sigma[type]*sigma[type]) ebond += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
    }

    // apply force to each of 2 atoms

    if (newton_bond || i1 < nlocal) {
      f[i1][0] += delx*fbond;
      f[i1][1] += dely*fbond;
      f[i1][2] += delz*fbond;
    }

    if (newton_bond || i2 < nlocal) {
      f[i2][0] -= delx*fbond;
      f[i2][1] -= dely*fbond;
      f[i2][2] -= delz*fbond;
    }

    if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);

    // subtract out pairwise contribution from 2 atoms via pair->single()
    // required since special_bond = 1,1,1
    // tally energy/virial in pair, using newton_bond as newton flag

    itype = atom->type[i1];
    jtype = atom->type[i2];

    if (rsq < cutsq[itype][jtype]) {
      evdwl = -force->pair->single(i1,i2,itype,jtype,rsq,1.0,1.0,fpair);
      fpair = -fpair;

      if (newton_bond || i1 < nlocal) {
	f[i1][0] += delx*fpair;
	f[i1][1] += dely*fpair;
	f[i1][2] += delz*fpair;
      }
      if (newton_bond || i2 < nlocal) {
	f[i2][0] -= delx*fpair;
	f[i2][1] -= dely*fpair;
	f[i2][2] -= delz*fpair;
      }

      if (evflag) force->pair->ev_tally(i1,i2,nlocal,newton_bond,
					evdwl,0.0,fpair,delx,dely,delz);
    }
  }
}