# Install/unInstall package classes in LAMMPS

if (test $1 = 1) then

    cp -p pair_reax_c.cpp ..

    cp -p fix_qeq_reax.cpp ..

    cp -p fix_reax_c.cpp ..

    cp -p pair_reax_c.h ..

    cp -p fix_qeq_reax.h ..

    cp -p fix_reax_c.h ..

    cp -p reaxc_allocate.cpp ..

    cp -p reaxc_basic_comm.cpp ..

    cp -p reaxc_bond_orders.cpp ..

    cp -p reaxc_bonds.cpp ..

    cp -p reaxc_control.cpp ..

    cp -p reaxc_ffield.cpp ..

    cp -p reaxc_forces.cpp ..

    cp -p reaxc_hydrogen_bonds.cpp ..

    cp -p reaxc_init_md.cpp ..

    cp -p reaxc_io_tools.cpp ..

    cp -p reaxc_list.cpp ..

    cp -p reaxc_lookup.cpp ..

    cp -p reaxc_multi_body.cpp ..

    cp -p reaxc_nonbonded.cpp ..

    cp -p reaxc_reset_tools.cpp ..

    cp -p reaxc_system_props.cpp ..

    cp -p reaxc_tool_box.cpp ..

    cp -p reaxc_torsion_angles.cpp ..

    cp -p reaxc_traj.cpp ..

    cp -p reaxc_valence_angles.cpp ..

    cp -p reaxc_vector.cpp ..

    cp -p reaxc_allocate.h ..

    cp -p reaxc_basic_comm.h ..

    cp -p reaxc_bond_orders.h ..

    cp -p reaxc_bonds.h ..

    cp -p reaxc_control.h ..

    cp -p reaxc_defs.h ..

    cp -p reaxc_ffield.h ..

    cp -p reaxc_forces.h ..

    cp -p reaxc_hydrogen_bonds.h ..

    cp -p reaxc_init_md.h ..

    cp -p reaxc_io_tools.h ..

    cp -p reaxc_list.h ..

    cp -p reaxc_lookup.h ..

    cp -p reaxc_multi_body.h ..

    cp -p reaxc_nonbonded.h ..

    cp -p reaxc_reset_tools.h ..

    cp -p reaxc_system_props.h ..

    cp -p reaxc_tool_box.h ..

    cp -p reaxc_torsion_angles.h ..

    cp -p reaxc_traj.h ..

    cp -p reaxc_types.h ..

    cp -p reaxc_valence_angles.h ..

    cp -p reaxc_vector.h ..

elif (test $1 = 0) then

    rm ../pair_reax_c.cpp

    rm ../fix_qeq_reax.cpp

    rm ../fix_reax_c.cpp

    rm ../pair_reax_c.h

    rm ../fix_qeq_reax.h

    rm ../fix_reax_c.h

    rm ../reaxc_allocate.cpp

    rm ../reaxc_basic_comm.cpp

    rm ../reaxc_bond_orders.cpp

    rm ../reaxc_bonds.cpp

    rm ../reaxc_control.cpp

    rm ../reaxc_ffield.cpp

    rm ../reaxc_forces.cpp

    rm ../reaxc_hydrogen_bonds.cpp

    rm ../reaxc_init_md.cpp

    rm ../reaxc_io_tools.cpp

    rm ../reaxc_list.cpp

    rm ../reaxc_lookup.cpp

    rm ../reaxc_multi_body.cpp

    rm ../reaxc_nonbonded.cpp

    rm ../reaxc_reset_tools.cpp

    rm ../reaxc_system_props.cpp

    rm ../reaxc_tool_box.cpp

    rm ../reaxc_torsion_angles.cpp

    rm ../reaxc_traj.cpp

    rm ../reaxc_valence_angles.cpp

    rm ../reaxc_vector.cpp

    rm ../reaxc_allocate.h

    rm ../reaxc_basic_comm.h

    rm ../reaxc_bond_orders.h

    rm ../reaxc_bonds.h

    rm ../reaxc_control.h

    rm ../reaxc_defs.h

    rm ../reaxc_ffield.h

    rm ../reaxc_forces.h

    rm ../reaxc_hydrogen_bonds.h

    rm ../reaxc_init_md.h

    rm ../reaxc_io_tools.h

    rm ../reaxc_list.h

    rm ../reaxc_lookup.h

    rm ../reaxc_multi_body.h

    rm ../reaxc_nonbonded.h

    rm ../reaxc_reset_tools.h

    rm ../reaxc_system_props.h

    rm ../reaxc_tool_box.h

    rm ../reaxc_torsion_angles.h

    rm ../reaxc_traj.h

    rm ../reaxc_types.h

    rm ../reaxc_valence_angles.h

    rm ../reaxc_vector.h

fi

The files in this directory are a user-contributed package for LAMMPS.

The person who created this package is Hasan Metin Aktulga, haktulga

at cs.purdue.edu, while at Purdue University.  Contact him directly,

or Aidan Thompson (Sandia) at athomps at sandia.gov, if you have

questions.

--------------------------------------

Note that the files with names starting with "reaxc" in this package

are from PuReMD, the Purdue ReaxFF Molecular Dynamics Program.  Its

copyright info and authorship info are listed below.

PACKAGE DESCRIPTION:

Contains a implementation for LAMMPS of the ReaxFF force field.

ReaxFF uses distance-dependent bond-order functions to represent the

contributions of chemical bonding to the potential energy.  It was

originally developed by Adri van Duin and the Goddard group at

CalTech.

The USER-REAXC version of ReaxFF (pair_style reax/c), implemented in

C++, should give identical or very similar results to pair_style reax,

which is a ReaxFF implementation on top of a Fortran library, a

version of which was originally authored by Adri van Duin.

The reax/c version should be somewhat faster and more scalable,

particularly with respect to the charge equilibration calculation.  It

should also be easier to build and use since there are no complicating

issues due to linking to a Fortran library.

OTHERS FILES INCLUDED:

User examples for pair_style reax/c are in examples/reax.

Thanks to Steve Plimpton and Aidan Thompson for their input on the

LAMMPS architecture and for their help in understanding and

customizing some of the required LAMMPS interfaces.

--------------------------------------

The reaxc files in this directory have the following header:

  PuReMD - Purdue ReaxFF Molecular Dynamics Program

  Copyright (2010) Purdue University

  Hasan Metin Aktulga, haktulga@cs.purdue.edu

  Joseph Fogarty, jcfogart@mail.usf.edu

  Sagar Pandit, pandit@usf.edu

  Ananth Y Grama, ayg@cs.purdue.edu

  This program is free software; you can redistribute it and/or

  modify it under the terms of the GNU General Public License as

  published by the Free Software Foundation; either version 2 of 

  the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,

  but WITHOUT ANY WARRANTY; without even the implied warranty of

  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  

  See the GNU General Public License for more details:

  <http://www.gnu.org/licenses/>.

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

   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(qeq/reax,FixQEqReax)

#else

#ifndef LMP_FIX_QEQ_REAX_H

#define LMP_FIX_QEQ_REAX_H

#include "fix.h"

namespace LAMMPS_NS {

class FixQEqReax : public Fix {

 public:

  FixQEqReax(class LAMMPS *, int, char **);

  ~FixQEqReax();

  int setmask();

  void init();

  void init_list(int,class NeighList *);

  void init_storage();

  void setup_pre_force(int);

  void pre_force(int);

  void setup_pre_force_respa(int, int);

  void pre_force_respa(int, int, int);

  void min_setup_pre_force(int);

  void min_pre_force(int);

  int matvecs;

  double qeq_time;

 private:

  int nevery,reaxflag;

  int n, N, m_fill;

  int n_cap, nmax, m_cap;

  int pack_flag;

  int nlevels_respa;

  class NeighList *list;

  double swa, swb;      // lower/upper Taper cutoff radius

  double Tap[8];        // Taper function

  double tolerance;     // tolerance for the norm of the rel residual in CG

  double *chi,*eta,*gamma;  // qeq parameters

  double **shld;

  // fictitious charges

  double *s, *t;

  double **s_hist, **t_hist;	

  int nprev;

  typedef struct{

    int n, m;

    int *firstnbr;

    int *numnbrs;

    int *jlist;

    double *val;

  } sparse_matrix;

  sparse_matrix H;

  double *Hdia_inv;

  double *b_s, *b_t;

  double *b_prc, *b_prm;

  //CG storage

  double *p, *q, *r, *d;

  //GMRES storage

  //double *g,*y;	  

  //double **v;

  //double **h;

  //double *hc, *hs;

  void pertype_parameters(char*);

  void init_shielding();

  void init_taper();

  void allocate_storage();

  void deallocate_storage();

  void reallocate_storage();

  void allocate_matrix();

  void deallocate_matrix();

  void reallocate_matrix();

  void init_matvec();

  void init_H();

  void compute_H();

  double calculate_H(double,double);

  void calculate_Q();

  int CG(double*,double*);

  //int GMRES(double*,double*);

  void sparse_matvec(sparse_matrix*,double*,double*);

  int pack_comm(int, int *, double *, int, int *);

  void unpack_comm(int, int, double *);

  int pack_reverse_comm(int, int, double *);

  void unpack_reverse_comm(int, int *, double *);

  double memory_usage();

  void grow_arrays(int);

  void copy_arrays(int, int);

  int pack_exchange(int, double *);

  int unpack_exchange(int, double *);

  double parallel_norm( double*, int );

  double parallel_dot( double*, double*, int );

  double parallel_vector_acc( double*, int );

  double norm(double*,int);

  void vector_sum(double*,double,double*,double,double*,int);

  void vector_scale(double*,double,double*,int);

  double dot(double*,double*,int);

  void vector_add(double*, double, double*,int);

};

}

#endif

#endif

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

   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: Hasan Metin Aktulga, Purdue University

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

#include "fix_reax_c.h"

#include "atom.h"

#include "pair.h"

#include "comm.h"

#include "memory.h"

using namespace LAMMPS_NS;

#define MAX_REAX_BONDS      30

#define MIN_REAX_BONDS      15

#define MIN_REAX_HBONDS     25

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

FixReaxC::FixReaxC(LAMMPS *lmp,int narg, char **arg) : 

  Fix(lmp, narg, arg)

{

  // perform initial allocation of atom-based arrays

  // register with atom class

  num_bonds = NULL;

  num_hbonds = NULL;

  grow_arrays(atom->nmax);

  atom->add_callback(0);

  // initialize arrays to MIN so atom migration is OK the 1st time

  int nlocal = atom->nlocal;

  for (int i = 0; i < nlocal; i++)

    num_bonds[i] = num_hbonds[i] = MIN_REAX_BONDS;

  // set comm sizes needed by this fix

  comm_forward = 1;

}

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

FixReaxC::~FixReaxC()

{

  // unregister this fix so atom class doesn't invoke it any more

  atom->delete_callback(id,0);

  // delete locally stored arrays

  memory->sfree(num_bonds);

  memory->sfree(num_hbonds);

}

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

int FixReaxC::setmask()

{

  int mask = 0;

  return mask;

}

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

   memory usage of local atom-based arrays

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

double FixReaxC::memory_usage()

{

  int nmax = atom->nmax;

  double bytes = nmax * 2 * sizeof(int);

  return bytes;

}

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

   allocate local atom-based arrays

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

void FixReaxC::grow_arrays(int nmax)

{

  num_bonds = (int *) memory->srealloc(num_bonds,nmax*sizeof(int),

				       "reaxc:num_bonds");

  num_hbonds = (int *) memory->srealloc(num_hbonds,nmax*sizeof(int),

					"reaxc:num_hbonds");

}

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

   copy values within local atom-based arrays

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

void FixReaxC::copy_arrays(int i, int j)

{

  num_bonds[j] = num_bonds[i];

  num_hbonds[j] = num_hbonds[i];

}

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

   pack values in local atom-based arrays for exchange with another proc

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

int FixReaxC::pack_exchange(int i, double *buf)

{

  buf[0] = num_bonds[i];

  buf[1] = num_hbonds[i];

  return 2;

}

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

   unpack values in local atom-based arrays from exchange with another proc

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

int FixReaxC::unpack_exchange(int nlocal, double *buf)

{

  num_bonds[nlocal] = static_cast<int> (buf[0]);

  num_hbonds[nlocal] = static_cast<int> (buf[1]);

  return 2;

}

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

int FixReaxC::pack_comm(int n, int *list, double *buf,

			 int pbc_flag, int *pbc)

{

  int i,j,m;

  m = 0;

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

    j = list[i];

    buf[m++] = num_bonds[j];

  }

  return 1;

}

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

void FixReaxC::unpack_comm(int n, int first, double *buf)

{

  int i,m,last;

  m = 0;

  last = first + n;

  for (i = first; i < last; i++)

    num_bonds[i] = buf[m++];

}