/* ----------------------------------------------------------------------
   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.
------------------------------------------------------------------------- */

#include "math.h"
#include "mpi.h"
#include "string.h"
#include "stdlib.h"
#include "fix_bond_update.h"
#include "update.h"
#include "respa.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "domain.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;
using namespace FixConst;

#define DELTA 16

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

FixBondUpdate::FixBondUpdate(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg)
{
  if (narg != 6 ) error->all(FLERR,"Illegal fix bond/break command");

  MPI_Comm_rank(world,&me);
  MPI_Comm_size(world,&nprocs);

  nevery = force->inumeric(FLERR,arg[3]);
  if (nevery <= 0) error->all(FLERR,"Illegal fix bond/break command");

  force_reneighbor = 1;
  next_reneighbor = -1;
  vector_flag = 1;
  size_vector = 2;
  global_freq = 1;
  extvector = 0;

  btype = force->inumeric(FLERR,arg[4]);
  cutoff = force->numeric(FLERR,arg[5]);

  if (btype < 1 || btype > atom->nbondtypes)
    error->all(FLERR,"Invalid bond type in fix bond/break command");
  if (cutoff < 0.0) error->all(FLERR,"Illegal fix bond/break command");

  cutsq = cutoff*cutoff;

  // error check

  if (atom->molecular != 1)
    error->all(FLERR,"Cannot use fix bond/break with non-molecular systems");

}

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

FixBondUpdate::~FixBondUpdate()
{
  delete random;

  // delete locally stored arrays

  memory->destroy(partner);
  memory->destroy(broken);
  delete [] copy;
}

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

int FixBondUpdate::setmask()
{
  int mask = 0;
  mask |= POST_INTEGRATE;
  return mask;
}

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

void FixBondUpdate::init()
{
  lastcheck = -1;

}

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

void FixBondUpdate::post_integrate()
{
  int i,j,ji,k,m,n,i1,i2,n1,n3,type;
  double delx,dely,delz,rsq;
  tagint *slist;

  if (update->ntimestep % nevery) return;

  int nlocal = atom->nlocal;
  int nall = atom->nlocal + atom->nghost;

  double **x = atom->x;
  tagint *tag = atom->tag;
  int *mask = atom->mask;

  int **bond_type = atom->bond_type;
  tagint **bond_atom = atom->bond_atom;
  int *num_bond = atom->num_bond;

  printf ("\nbefore the main loop %d ",update->ntimestep  );
  for (i = 0; i < nlocal; i++) {
	
    // delete bond from atom I if I stores it
    // atom J will also do this
    //printf ("\nloop in atom # %d with tag %d", i,tag[i]);

    for (m = 0; m < num_bond[i]; m++) {
	    if (!(mask[i] & groupbit)) continue;
	    
	    if (abs(bond_type[i][m]) != btype) continue;
		
	    // j is the tag of partner atom
	    // ji is the current index of the partner atom
	    
	    j = bond_atom[i][m]; 
	    ji = atom->map(j);
	    
	    if (!(mask[ji] & groupbit)) continue;

	    // check if we have j in this domain
	    // it is ok if we don't, it means they are far away
	    if (ji == -1) continue; 

	    delx = x[i][0] - x[ji][0];
	    dely = x[i][1] - x[ji][1];
	    delz = x[i][2] - x[ji][2];

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

	    if (bond_type[i][m] == btype) {
		    // break the bond
		    if (rsq > cutsq)   
			    disable_bond(i,m);
	    }
	    
	    else if (bond_type[i][m] == -btype) {
		    // un-break the bond
		    if (rsq <= cutsq)
			    enable_bond (i,m);
	    }
	
        }
  }

//  next_reneighbor = update->ntimestep;
}

/*------------
disable the Bth bond and all dihedrals owned by atom M
*/

void FixBondUpdate::disable_bond(int m, int b) {
	
	int i;
	tagint partner;

	int *bond_type = atom->bond_type[m];
  	tagint *bond_atom = atom->bond_atom[m];
	
	int num_dihedral = atom->num_dihedral[m];
	int *dihedral_type = atom->dihedral_type[m];

	printf("\nstep : %d, disable atom %d bond %d with %d dihedrals.",update->ntimestep, m,b,num_dihedral);
 
	// disable the jth bond
	bond_type[b] = -btype;
	
	// disable all dihedrals
	for (i = 0; i < num_dihedral; i++) 
		dihedral_type[i] = -abs(dihedral_type[i]);

}

/*------------
enable the Bth bond and all dihedrals owned by atom M
*/

void FixBondUpdate::enable_bond(int m, int b) {
	printf("enabling atom %d bond %d", m,b);
	
	int i;
	tagint partner;

	int *bond_type = atom->bond_type[m];
  	tagint *bond_atom = atom->bond_atom[m];
	
	int num_dihedral = atom->num_dihedral[m];
	int *dihedral_type = atom->dihedral_type[m];

	printf("\nstep : %d, enabling atom %d bond %d with %d dihedrals.",update->ntimestep, m,b,num_dihedral);
 
	// enable the jth bond
	bond_type[b] = btype;
	
	// enable all dihedrals
	for (i = 0; i < num_dihedral; i++) 
		dihedral_type[i] = abs(dihedral_type[i]);

}