/* -*- c++ -*- ----------------------------------------------------------
   This command computes the number of collisions per atom based on a 
   separation distance. This distance is computed as the sum of the 
   individual radiuses of the atoms. The ReaxFF pair style is required.

   Collisions are detected as follows:
   a) If an atom has the same molecule ID with the considered atom (i.e. 
      an atom that belong on the same molecule), no collision is counted
   b) A collision with an atom is broadcasted in all atoms of the same 
      molecule (i.e. between a collision pair all atoms acquire that 
      information)
   c) If an atom in the neighbor list has a molecule ID that was already
      considered as collision, no collision is counted again

   At the end all atoms of the same molecule have the total number of collisions
   for that particular molecule, i.e. if a C12H26 molecule collides with 3 O2 molecules
   all atoms of the C12H26 molecule will have a count of collisions equal to 3
   
   Special treatment for the ghost atoms of a molecule that is located between two processors
   was added. The number of collisions is summed over all atoms of the molecule, even if that
   is located in different processors.

   NOTE: If the reax/c/species consinders an interaction close enough that a bond is 
   created, then this may not be counted as a collision, since all atoms of those
   involved molecules will have the same molecular ID

   Structure:
    
      compute ID group-ID collision/contact type1 radius1 type2 radius2 ... flag value
   
   - ID, group-ID are documented in compute command
   - temp/comp/chunk = style name of this compute command
   - value = type1 radius1 type2 radius2 ...
           - type1   : type of atom 1 (in integer format)
           - radius1 : radius of collision of type 1 
           - type2   : type of atom 2 (in integer format)
           - radius2 : radius of collision of type 2 
   - flag = start
           - start : start timestep that the chunk takes values,
                     ex. if chunk based on fix reax/c/species 1 1 100 then
                     before dt=100 the compute chunk will have no values thus
                     this compute will show arbitrary values or give a segmentation
                     memory fault
   - value = starting timestep 


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

#ifdef COMPUTE_CLASS

ComputeStyle(collision/contact,ComputeCollisionContact)

#else

#ifndef LMP_COMPUTE_COLLISION_CONTACT_H
#define LMP_COMPUTE_COLLISION_CONTACT_H

#include "compute.h"
#include "fix.h"
#include "reaxc_types.h"
#include "pair_reaxc.h"
#include "reaxc_defs.h"
#include "fix_reaxc_species.h"

namespace LAMMPS_NS {

class ComputeCollisionContact : public Compute {
 public:
  ComputeCollisionContact(class LAMMPS *, int, char **);
  ~ComputeCollisionContact();
  void init();
  void init_list(int, class NeighList *);
  void compute_peratom();
  int pack_reverse_comm(int, int, double *);
  void unpack_reverse_comm(int, int *, double *);
  double memory_usage();

 private:
  int nmax, nchunk, nmol;
  int startstep;
  class NeighList *list;
  class PairReaxC *reaxc;
  double *contact, *contactall, *contactPerAtom;
  double *radiuses;
  int *atomChunkID;
  int **neighbourChunk;
  int me,nprocs;
};

}

#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: Compute collision/contact requires a pair style be defined

Self-explanatory.

W: More than one compute collision/contact

It is not efficient to use compute collision/contact more than once.

*/