/* ----------------------------------------------------------------------
   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 "mpi.h"
#include "stdlib.h"
#include "string.h"
#include "compute_SW_test.h"
#include "atom.h"
#include "force.h"
#include "modify.h"
#include "comm.h"
#include "memory.h"
#include "group.h"
#include "error.h"
#include "update.h"
#include "neighbor.h"
#include "math.h"
#include "pair.h"
#include "universe.h"

using namespace LAMMPS_NS;

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

ComputeSWTest::ComputeSWTest(LAMMPS *lmp, int narg, char **arg) : 
  Compute(lmp, narg, arg)
{
  if (narg != 3) error->all("Illegal compute SW/test command");

  scalar_flag = 1;
  vector_flag = 1;
  size_vector = 4;
  
  comm_reverse = 1;
  neigh_half_once = 1;

  extensive = 0;
  energy1 = NULL;
  nmax =0;
  vector = new double[4];

}

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

ComputeSWTest::~ComputeSWTest()
{ 
  memory->sfree(energy1);
  delete [] vector;
}

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

void ComputeSWTest::init()
{

}

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

void ComputeSWTest::compute_vector()
{
  int i,j,k,m,n,itag,jtag,itype,jtype,ktype,iparam,numneigh;
  double xtmp,ytmp,ztmp,delx,dely,delz;
  double rsq,rsq1,rsq2,eng,fforce;
  double delr1[3],delr2[3],fj[3],fk[3];
  int *neighs;
  double **f;


  // allocate memory for the vector energy1
  if (atom->nmax > nmax) {
    memory->sfree(energy1);
    nmax = atom->nmax;
    energy1 = (double *)
      memory->smalloc(nmax*sizeof(double),"compute/epair/atom:energy");
    scalar_atom = energy1;
  }

  // set elements of energy1 to zero
  if (force->newton_pair) n = atom->nlocal + atom->nghost;
  else n = atom->nlocal;
  for (i = 0; i < n; i++) energy1[i] = 0.0;
  
  if (!neighbor->half_every) neighbor->build_half();
  
  double **cutsq = force->pair->cutsq;

  double **x = atom->x;
  double **v = atom->v; 

  int *tag = atom->tag;
  int *type = atom->type;
  int nlocal = atom->nlocal;

  Pair::Param params;  
  double t[4];


  // loop over full neighbor list of my atoms
  
  for (i = 0; i < nlocal; i++) {
    itag = tag[i];
    // next line leads to memory access failure
    itype = force->pair->map[type[i]];
    xtmp = x[i][0];
    ytmp = x[i][1];
    ztmp = x[i][2];

    // two-body interactions

    neighs = neighbor->firstneigh_full[i];
    numneigh = neighbor->numneigh_full[i];

    for (m = 0; m < numneigh; m++) {

      j = neighs[m];
      jtag = tag[j];
      jtype = force->pair->map[type[j]];
      
      delx = xtmp - x[j][0];
      dely = ytmp - x[j][1];
      delz = ztmp - x[j][2];
      rsq = delx*delx + dely*dely + delz*delz;
      iparam = force->pair->elem2param[itype][jtype][jtype];
      if (rsq > force->pair->params[iparam].cutsq) continue;
      force->pair->twobody(&force->pair->params[iparam],rsq,fforce,1,eng); 
      %%% Diagnostic code deleted for simplicity

    } // end of loop over m (two-body interactions)

    // three-body interactions
    for (m = 0; m < numneigh-1; m++) {
      j = neighs[m];
      jtype = force->pair->map[type[j]];

      for (n = m+1; n < numneigh; n++) {
	k = neighs[n];
	ktype = force->pair->map[type[k]];
	iparam = force->pair->elem2param[itype][jtype][ktype];

	delr1[0] = x[j][0] - xtmp;
	delr1[1] = x[j][1] - ytmp;
	delr1[2] = x[j][2] - ztmp;
	rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
	if (rsq1 > force->pair->params[iparam].cutsq) continue;

	delr2[0] = x[k][0] - xtmp;
	delr2[1] = x[k][1] - ytmp;
	delr2[2] = x[k][2] - ztmp;
	rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
	if (rsq2 > force->pair->params[iparam].cutsq) continue;
	force->pair->threebody(&force->pair->params[iparam],rsq1,rsq2,delr1,delr2,fj,fk,1,eng);
        %%% Diagnostic code deleted for simplicity

      } // end of loop over n
    } // end of loop over m (3-body interactions)
    
  } // end of loop over i (my atoms)
 
  %%% Diagnostic code deleted for simplicity

  
  MPI_Allreduce(t,vector,4,MPI_DOUBLE,MPI_SUM,world);
}


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

int ComputeSWTest::pack_reverse_comm(int n, int first, double *buf)
{
  int i,m,last;
  
  m = 0;
  last = first + n;
  for (i = first; i < last; i++) buf[m++] = energy1[i];
  return 1;
}

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

void ComputeSWTest::unpack_reverse_comm(int n, int *list, double *buf)
{
  int i,j,m;
  m = 0;
  for (i = 0; i < n; i++) {
    j = list[i];
    energy1[j] += buf[m++];
  }
}