/* ----------------------------------------------------------------------
   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: Alexander Stukowski
                        Technical University of Darmstadt,
                        Germany Department of Materials Science
------------------------------------------------------------------------- */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_cdeam.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

// This is for debugging purposes. The ASSERT() macro is used in the code to check
// if everything runs as expected. Change this to #if 0 if you don't need the checking.
#if 0
        #define ASSERT(cond) ((!(cond)) ? my_failure(error,__FILE__,__LINE__) : my_noop())

        inline void my_noop() {}
        inline void my_failure(Error* error, const char* file, int line) {
                char str[1024];
                sprintf(str,"Assertion failure: File %s, line %i", file, line);
                error->one(FLERR,str);
        }
#else
        #define ASSERT(cond)
#endif

#define MAXLINE 1024        // This sets the maximum line length in EAM input files.

PairCDEAM::PairCDEAM(LAMMPS *lmp, int _cdeamVersion) : PairEAM(lmp), PairEAMAlloy(lmp), cdeamVersion(_cdeamVersion)
{
        single_enable = 0;
        restartinfo = 0;

        rhoB = NULL;
        D_values = NULL;
        hcoeff = NULL;

        // Set communication buffer sizes needed by this pair style.
        if(cdeamVersion == 1) {
                comm_forward = 4;
                comm_reverse = 3;
        }
        else if(cdeamVersion == 2) {
                comm_forward = 3;
                comm_reverse = 2;
        }
        else {
                error->all(FLERR,"Invalid CD-EAM potential version.");
        }
}

PairCDEAM::~PairCDEAM()
{
        memory->destroy(rhoB);
        memory->destroy(D_values);
        if(hcoeff) delete[] hcoeff;
}

void PairCDEAM::compute(int eflag, int vflag)
{
        int i,j,ii,jj,inum,jnum,itype,jtype;
        double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
        double rsq,rhoip,rhojp,recip,phi;
        int *ilist,*jlist,*numneigh,**firstneigh;

        evdwl = 0.0;
        if (eflag || vflag) ev_setup(eflag,vflag);
        else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;

        // Grow per-atom arrays if necessary
        if(atom->nmax > nmax) {
                memory->destroy(rho);
                memory->destroy(fp);
                memory->destroy(rhoB);
                memory->destroy(D_values);
                nmax = atom->nmax;
                memory->create(rho,nmax,"pair:rho");
                memory->create(rhoB,nmax,"pair:rhoB");
                memory->create(fp,nmax,"pair:fp");
                memory->create(D_values,nmax,"pair:D_values");
        }

        double **x = atom->x;
        double **f = atom->f;
        int *type = atom->type;
        int nlocal = atom->nlocal;
        int newton_pair = force->newton_pair;

        inum = list->inum;
        ilist = list->ilist;
        numneigh = list->numneigh;
        firstneigh = list->firstneigh;

        // Zero out per-atom arrays.
        int m = nlocal + atom->nghost;
        for(i = 0; i < m; i++) {
                rho[i] = 0.0;
                rhoB[i] = 0.0;
                D_values[i] = 0.0;
        }

        // Stage I

        // Compute rho and rhoB at each local atom site.
        // Additionally calculate the D_i values here if we are using the one-site formulation.
        // For the two-site formulation we have to calculate the D values in an extra loop (Stage II).
        for(ii = 0; ii < inum; ii++) {
                i = ilist[ii];
                xtmp = x[i][0];
                ytmp = x[i][1];
                ztmp = x[i][2];
                itype = type[i];
                jlist = firstneigh[i];
                jnum = numneigh[i];

                for(jj = 0; jj < jnum; jj++) {
                        j = jlist[jj];
                        j &= NEIGHMASK;

                        delx = xtmp - x[j][0];
                        dely = ytmp - x[j][1];
                        delz = ztmp - x[j][2];
                        rsq = delx*delx + dely*dely + delz*delz;

                        if(rsq < cutforcesq) {
                                jtype = type[j];
                                double r = sqrt(rsq);
                                const EAMTableIndex index = radiusToTableIndex(r);
                                double localrho = RhoOfR(index, jtype, itype);
                                rho[i] += localrho;
                                if(jtype == speciesB) rhoB[i] += localrho;
                                if(newton_pair || j < nlocal) {
                                        localrho = RhoOfR(index, itype, jtype);
                                        rho[j] += localrho;
                                        if(itype == speciesB) rhoB[j] += localrho;
                                }

                                if(cdeamVersion == 1 && itype != jtype) {
                                        // Note: if the i-j interaction is not concentration dependent (because either
                                        // i or j are not species A or B) then its contribution to D_i and D_j should
                                        // be ignored.
                                        // This if-clause is only required for a ternary.
                                        if((itype == speciesA && jtype == speciesB) || (jtype == speciesA && itype == speciesB)) {
                                                double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);
                                                D_values[i] += Phi_AB;
                                                if(newton_pair || j < nlocal)
                                                        D_values[j] += Phi_AB;
                                        }
                                }
                        }
                }
        }

        // Communicate and sum densities.
        if(newton_pair) {
                communicationStage = 1;
                comm->reverse_comm_pair(this);
        }

        // fp = derivative of embedding energy at each atom
        // phi = embedding energy at each atom
        for(ii = 0; ii < inum; ii++) {
                i = ilist[ii];
                EAMTableIndex index = rhoToTableIndex(rho[i]);
                fp[i] = FPrimeOfRho(index, type[i]);
                if(eflag) {
                        phi = FofRho(index, type[i]);
                        if (eflag_global) eng_vdwl += phi;
                        if (eflag_atom) eatom[i] += phi;
                }
        }

        // Communicate derivative of embedding function and densities
        // and D_values (this for one-site formulation only).
        communicationStage = 2;
        comm->forward_comm_pair(this);

        // The electron densities may not drop to zero because then the concentration would no longer be defined.
        // But the concentration is not needed anyway if there is no interaction with another atom, which is the case
        // if the electron density is exactly zero. That's why the following lines have been commented out.
        //
        //for(i = 0; i < nlocal + atom->nghost; i++) {
        //        if(rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB))
        //                error->one(FLERR,"CD-EAM potential routine: Detected atom with zero electron density.");
        //}

        // Stage II
        // This is only required for the original two-site formulation of the CD-EAM potential.

        if(cdeamVersion == 2) {
                // Compute intermediate value D_i for each atom.
                for(ii = 0; ii < inum; ii++) {
                        i = ilist[ii];
                        xtmp = x[i][0];
                        ytmp = x[i][1];
                        ztmp = x[i][2];
                        itype = type[i];
                        jlist = firstneigh[i];
                        jnum = numneigh[i];

                        // This code line is required for ternary alloys.
                        if(itype != speciesA && itype != speciesB) continue;

                        double x_i = rhoB[i] / rho[i];        // Concentration at atom i.

                        for(jj = 0; jj < jnum; jj++) {
                                j = jlist[jj];
                                j &= NEIGHMASK;
                                jtype = type[j];
                                if(itype == jtype) continue;

                                // This code line is required for ternary alloys.
                                if(jtype != speciesA && jtype != speciesB) continue;

                                delx = xtmp - x[j][0];
                                dely = ytmp - x[j][1];
                                delz = ztmp - x[j][2];
                                rsq = delx*delx + dely*dely + delz*delz;

                                if(rsq < cutforcesq) {
                                        double r = sqrt(rsq);
                                        const EAMTableIndex index = radiusToTableIndex(r);

                                        // The concentration independent part of the cross pair potential.
                                        double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r);

                                        // Average concentration of two sites
                                        double x_ij = 0.5 * (x_i + rhoB[j]/rho[j]);

                                        // Calculate derivative of h(x_ij) polynomial function.
                                        double h_prime = evalHprime(x_ij);

                                        D_values[i] += h_prime * Phi_AB / (2.0 * rho[i] * rho[i]);
                                        if(newton_pair || j < nlocal)
                                                D_values[j] += h_prime * Phi_AB / (2.0 * rho[j] * rho[j]);
                                }
                        }
                }

                // Communicate and sum D values.
                if(newton_pair) {
                        communicationStage = 3;
                        comm->reverse_comm_pair(this);
                }
                communicationStage = 4;
                comm->forward_comm_pair(this);
        }

        // Stage III

        // Compute force acting on each atom.
        for(ii = 0; ii < inum; ii++) {
                i = ilist[ii];
                xtmp = x[i][0];
                ytmp = x[i][1];
                ztmp = x[i][2];
                itype = type[i];

                jlist = firstneigh[i];
                jnum = numneigh[i];

                // Concentration at site i
                double x_i = -1.0;                // The value -1 indicates: no concentration dependence for all interactions of atom i.
                                                                // It will be replaced by the concentration at site i if atom i is either A or B.

                double D_i, h_prime_i;

                // This if-clause is only required for ternary alloys.
                if((itype == speciesA || itype == speciesB) && rho[i] != 0.0) {

                        // Compute local concentration at site i.
                        x_i = rhoB[i]/rho[i];
                        ASSERT(x_i >= 0 && x_i<=1.0);

                        if(cdeamVersion == 1) {
                                // Calculate derivative of h(x_i) polynomial function.
                                h_prime_i = evalHprime(x_i);
                                D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]);
                        } else if(cdeamVersion == 2) {
                                D_i = D_values[i];
                        } else {
                          ASSERT(false);
                        }
                }

                for(jj = 0; jj < jnum; jj++) {
                        j = jlist[jj];
                        j &= NEIGHMASK;

                        delx = xtmp - x[j][0];
                        dely = ytmp - x[j][1];
                        delz = ztmp - x[j][2];
                        rsq = delx*delx + dely*dely + delz*delz;

                        if(rsq < cutforcesq) {
                                jtype = type[j];
                                double r = sqrt(rsq);
                                const EAMTableIndex index = radiusToTableIndex(r);

                                // rhoip = derivative of (density at atom j due to atom i)
                                // rhojp = derivative of (density at atom i due to atom j)
                                // psip needs both fp[i] and fp[j] terms since r_ij appears in two
                                //   terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
                                //   hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
                                rhoip = RhoPrimeOfR(index, itype, jtype);
                                rhojp = RhoPrimeOfR(index, jtype, itype);
                                fpair = fp[i]*rhojp + fp[j]*rhoip;
                                recip = 1.0/r;

                                double x_j = -1;  // The value -1 indicates: no concentration dependence for this i-j pair
                                                  // because atom j is not of species A nor B.

                                // This code line is required for ternary alloy.
                                if(jtype == speciesA || jtype == speciesB) {
                                        ASSERT(rho[i] != 0.0);
                                        ASSERT(rho[j] != 0.0);

                                        // Compute local concentration at site j.
                                        x_j = rhoB[j]/rho[j];
                                        ASSERT(x_j >= 0 && x_j<=1.0);

                                        double D_j=0.0;
                                        if(cdeamVersion == 1) {
                                                // Calculate derivative of h(x_j) polynomial function.
                                                double h_prime_j = evalHprime(x_j);
                                                D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]);
                                        } else if(cdeamVersion == 2) {
                                                D_j = D_values[j];
                                        } else {
                                          ASSERT(false);
                                        }
                                        double t2 = -rhoB[j];
                                        if(itype == speciesB) t2 += rho[j];
                                        fpair += D_j * rhoip * t2;
                                }

                                // This if-clause is only required for a ternary alloy.
                                // Actually we don't need it at all because D_i should be zero anyway if
                                // atom i has no concentration dependent interactions (because it is not species A or B).
                                if(x_i != -1.0) {
                                        double t1 = -rhoB[i];
                                        if(jtype == speciesB) t1 += rho[i];
                                        fpair += D_i * rhojp * t1;
                                }

                                double phip;
                                double phi = PhiOfR(index, itype, jtype, recip, phip);
                                if(itype == jtype || x_i == -1.0 || x_j == -1.0) {
                                        // Case of no concentration dependence.
                                        fpair += phip;
                                } else {
                                        // We have a concentration dependence for the i-j interaction.
                                        double h=0.0;
                                        if(cdeamVersion == 1) {
                                                // Calculate h(x_i) polynomial function.
                                                double h_i = evalH(x_i);
                                                // Calculate h(x_j) polynomial function.
                                                double h_j = evalH(x_j);
                                                h = 0.5 * (h_i + h_j);
                                        } else if(cdeamVersion == 2) {
                                                // Average concentration.
                                                double x_ij = 0.5 * (x_i + x_j);
                                                // Calculate h(x_ij) polynomial function.
                                                h = evalH(x_ij);
                                        } else {
                                          ASSERT(false);
                                        }
                                        fpair += h * phip;
                                        phi *= h;
                                }

                                // Divide by r_ij and negate to get forces from gradient.
                                fpair /= -r;

                                f[i][0] += delx*fpair;
                                f[i][1] += dely*fpair;
                                f[i][2] += delz*fpair;
                                if(newton_pair || j < nlocal) {
                                        f[j][0] -= delx*fpair;
                                        f[j][1] -= dely*fpair;
                                        f[j][2] -= delz*fpair;
                                }

                                if(eflag) evdwl = phi;
                                if(evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
                        }
                }
        }

        if(vflag_fdotr) virial_fdotr_compute();
}

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

void PairCDEAM::coeff(int narg, char **arg)
{
        PairEAMAlloy::coeff(narg, arg);

        // Make sure the EAM file is a CD-EAM binary alloy.
        if(setfl->nelements < 2)
                error->all(FLERR,"The EAM file must contain at least 2 elements to be used with the eam/cd pair style.");

        // Read in the coefficients of the h polynomial from the end of the EAM file.
        read_h_coeff(arg[2]);

        // Determine which atom type is the A species and which is the B species in the alloy.
        // By default take the first element (index 0) in the EAM file as the A species
        // and the second element (index 1) in the EAM file as the B species.
        speciesA = -1;
        speciesB = -1;
        for(int i = 1; i <= atom->ntypes; i++) {
                if(map[i] == 0) {
                        if(speciesA >= 0)
                                error->all(FLERR,"The first element from the EAM file may only be mapped to a single atom type.");
                        speciesA = i;
                }
                if(map[i] == 1) {
                        if(speciesB >= 0)
                                error->all(FLERR,"The second element from the EAM file may only be mapped to a single atom type.");
                        speciesB = i;
                }
        }
        if(speciesA < 0)
                error->all(FLERR,"The first element from the EAM file must be mapped to exactly one atom type.");
        if(speciesB < 0)
                error->all(FLERR,"The second element from the EAM file must be mapped to exactly one atom type.");
}

/* ----------------------------------------------------------------------
   Reads in the h(x) polynomial coefficients
------------------------------------------------------------------------- */
void PairCDEAM::read_h_coeff(char *filename)
{
        if(comm->me == 0) {
                // Open potential file
                FILE *fptr;
                char line[MAXLINE];
                char nextline[MAXLINE];
                fptr = force->open_potential(filename);
                if (fptr == NULL) {
                        char str[128];
                        sprintf(str,"Cannot open EAM potential file %s", filename);
                        error->one(FLERR,str);
                }

                // h coefficients are stored at the end of the file.
                // Skip to last line of file.
                while(fgets(nextline, MAXLINE, fptr) != NULL) {
                        strcpy(line, nextline);
                }
                char* ptr = strtok(line, " \t\n\r\f");
                int degree = atoi(ptr);
                nhcoeff = degree+1;
                hcoeff = new double[nhcoeff];
                int i = 0;
                while((ptr = strtok(NULL," \t\n\r\f")) != NULL && i < nhcoeff) {
                        hcoeff[i++] = atof(ptr);
                }
                if(i != nhcoeff || nhcoeff < 1)
                        error->one(FLERR,"Failed to read h(x) function coefficients from EAM file.");

                // Close the potential file.
                fclose(fptr);
        }

        MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world);
        if(comm->me != 0) hcoeff = new double[nhcoeff];
        MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world);
}


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

int PairCDEAM::pack_forward_comm(int n, int *list, double *buf,
                                 int pbc_flag, int *pbc)
{
        int i,j,m;

        m = 0;
        if(communicationStage == 2) {
                if(cdeamVersion == 1) {
                        for (i = 0; i < n; i++) {
                                j = list[i];
                                buf[m++] = fp[j];
                                buf[m++] = rho[j];
                                buf[m++] = rhoB[j];
                                buf[m++] = D_values[j];
                        }
                        return m;
                }
                else if(cdeamVersion == 2) {
                        for (i = 0; i < n; i++) {
                                j = list[i];
                                buf[m++] = fp[j];
                                buf[m++] = rho[j];
                                buf[m++] = rhoB[j];
                        }
                        return m;
                }
                else { ASSERT(false); return 0; }
        }
        else if(communicationStage == 4) {
                for (i = 0; i < n; i++) {
                        j = list[i];
                        buf[m++] = D_values[j];
                }
                return m;
        }
        else return 0;
}

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

void PairCDEAM::unpack_forward_comm(int n, int first, double *buf)
{
        int i,m,last;

        m = 0;
        last = first + n;
        if(communicationStage == 2) {
                if(cdeamVersion == 1) {
                        for(i = first; i < last; i++) {
                                fp[i] = buf[m++];
                                rho[i] = buf[m++];
                                rhoB[i] = buf[m++];
                                D_values[i] = buf[m++];
                        }
                }
                else if(cdeamVersion == 2) {
                        for(i = first; i < last; i++) {
                                fp[i] = buf[m++];
                                rho[i] = buf[m++];
                                rhoB[i] = buf[m++];
                        }
                } else {
                  ASSERT(false);
                }
        }
        else if(communicationStage == 4) {
                for(i = first; i < last; i++) {
                        D_values[i] = buf[m++];
                }
        }
}

/* ---------------------------------------------------------------------- */
int PairCDEAM::pack_reverse_comm(int n, int first, double *buf)
{
        int i,m,last;

        m = 0;
        last = first + n;

        if(communicationStage == 1) {
                if(cdeamVersion == 1) {
                        for(i = first; i < last; i++) {
                                buf[m++] = rho[i];
                                buf[m++] = rhoB[i];
                                buf[m++] = D_values[i];
                        }
                        return m;
                }
                else if(cdeamVersion == 2) {
                        for(i = first; i < last; i++) {
                                buf[m++] = rho[i];
                                buf[m++] = rhoB[i];
                        }
                        return m;
                }
                else { ASSERT(false); return 0; }
        }
        else if(communicationStage == 3) {
                for(i = first; i < last; i++) {
                        buf[m++] = D_values[i];
                }
                return m;
        }
        else return 0;
}

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

void PairCDEAM::unpack_reverse_comm(int n, int *list, double *buf)
{
        int i,j,m;

        m = 0;
        if(communicationStage == 1) {
                if(cdeamVersion == 1) {
                        for(i = 0; i < n; i++) {
                                j = list[i];
                                rho[j] += buf[m++];
                                rhoB[j] += buf[m++];
                                D_values[j] += buf[m++];
                        }
                } else if(cdeamVersion == 2) {
                        for(i = 0; i < n; i++) {
                                j = list[i];
                                rho[j] += buf[m++];
                                rhoB[j] += buf[m++];
                        }
                } else {
                  ASSERT(false);
                }
        }
        else if(communicationStage == 3) {
                for(i = 0; i < n; i++) {
                        j = list[i];
                        D_values[j] += buf[m++];
                }
        }
}

/* ----------------------------------------------------------------------
   memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairCDEAM::memory_usage()
{
        double bytes = 2 * nmax * sizeof(double);
        return PairEAMAlloy::memory_usage() + bytes;
}
