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

FixStyle(rattle,FixRattle)

#else

#ifndef LMP_FIX_RATTLE_H
#define LMP_FIX_RATTLE_H

#include "fix.h"
#include "fix_shake.h"
namespace LAMMPS_NS {


// For testing and debugging. Set RATTLE_TEST to 1 if you want the contraints to be checked
// at the end of the time step. Set RATTLE_DEBUG to 1 if you want some more output for debuggin purposes.
#define RATTLE_DEBUG 0
#define RATTLE_TEST  1



enum COMM_MODE {V, VP, XSHAKE};

class FixRattle : public FixShake {

 double ** vp;                      // for unconstrained velocity update 
 double dt, dtfvo2, dtfsq_full;     // timesteps with conversion factors
 COMM_MODE comm_mode;               // determines mode for forwad_packing and forward_unpacking

 int pid;                           // DEBUG: processor id
 FILE * dfp;                        // DEBUG: stream



 public:
  FixRattle(class LAMMPS *, int, char **);
  ~FixRattle();
  virtual void post_force(int);
  virtual void initial_integrate(int);
  void vrattle(int m);
  void vrattle3(int m);
  void vrattle4(int m);
  void vrattle3angle(int m);
  void solve3x3exactly(const double a[][3], const double c[], double l[]);
  void update_v_half_nocons();
  int  setmask();

  virtual double memory_usage();
  virtual void grow_arrays(int);
  
  virtual int pack_forward_comm(int n, int *list, double *buf,int pbc_flag, int *pbc);
  virtual void unpack_forward_comm(int n, int first, double *buf);


  // DEBUG
  bool check_molecule(double ** v, int m, bool checkr, bool checkv);
  bool check_constraints(double **v, bool checkr, bool checkv);
  void end_of_step();

};
}
#endif
#endif