refine the code and documentation

\documentclass[aps,pr,onecolumn,superscriptaddress,noshowpacs,a4paper,15pt]{revtex4}

\pdfoutput=1

\bibliographystyle{apsrev4}

\usepackage{color}

\usepackage{dcolumn} %Align table columns on decimal point

\usepackage{amssymb}

\usepackage{amsmath}

\usepackage{amsthm}

\usepackage{graphicx}

\usepackage[pdftex]{hyperref}

\hypersetup{colorlinks=true,citecolor=blue,linkcolor=red,urlcolor=blue}

\usepackage[all]{hypcap}

\newcommand{\red}{\color{red}}

\newcommand{\blue}{\color{blue}}

\definecolor{green}{rgb}{0,0.5,0}

\newcommand{\green}{\color{green}}

\newcommand{\white}{\color{white}}

%\newcommand{\cite}[1]{\hspace{-1 ex} % \nocite{#1}\citenum{#1}}

\thickmuskip=0.5\thickmuskip %shorter spaces in math

%

\begin{document}

%

\begingroup

\Large

\begin{eqnarray*}

  E & = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\[15pt]

  V_{ij} & = & {\rm Tap}(r_{ij})\left \{ e^{-\alpha (r_{ij}/\beta -1)} 

               \left [ \epsilon + f(\rho_{ij}) + f(\rho_{ji})\right ] - 

                \frac{1}{1+e^{-d\left [ \left ( r_{ij}/\left (s_R \cdot r^{eff} \right ) \right )-1 \right ]}}

                \cdot \frac{C_6}{r^6_{ij}} \right \}\\[15pt]

  \rho_{ij}^2 & = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_i)^2 \\[15pt]

  \rho_{ji}^2 & = & r_{ij}^2 - ({\bf r}_{ij} \cdot {\bf n}_j)^2 \\[15pt]

  f(\rho) & = &  C e^{ -( \rho / \delta )^2 }\\[15pt]

  {\rm Tap}(r_{ij}) & = & 20\left ( \frac{r_{ij}}{R_{cut}} \right )^7 -

                          70\left ( \frac{r_{ij}}{R_{cut}} \right )^6 +

                          84\left ( \frac{r_{ij}}{R_{cut}} \right )^5 -

                          35\left ( \frac{r_{ij}}{R_{cut}} \right )^4 + 1

\end{eqnarray*}

\endgroup

%

\end{document}

%

[Description:]

The {ILP/graphene/hBN} style computes the registry-dependent interlayer

The {ilp/graphene/hbn} style computes the registry-dependent interlayer

potential (RDILP) potential as described in "(Leven)"_#Leven and

"(Maaravi)"_#Maaravi2. The normals are calculated in the way as described

in "(Kolmogorov)"_#Kolmogorov2.

:c,image(Eqs/pair_ilp_gr_hBN.jpg)

:c,image(Eqs/pair_ilp_graphene_hbn.jpg)

Where Tap(r_ij) is the taper function which provides a continuous

cutoff (up to third derivative) for interatomic separations larger than

# Kolmogorov-Crespi Potential

#

# Cite as A.N. Kolmogorov & V. H. Crespi, 

# Registry-dependent interlayer potential for graphitic systems

# Physical Review B 71, 235415 (2005)

# 

#      z0     C0     C2      C4       C    delta  lambda      A     S	rcut

C  C  3.34  15.71  12.29   4.933   3.030   0.578   3.629   10.238  1.0  2.0

    memory->destroy(dnormal);

    memory->destroy(dnormdri);

    nmax = atom->nmax;

    memory->create(normal,nmax,3,"ILP:normal");

    memory->create(dnormdri,3,3,nmax,"ILP:dnormdri");

    memory->create(dnormal,3,3,3,nmax,"ILP:dnormal");

    memory->create(normal,nmax,3,"ILPGrapheneHBN:normal");

    memory->create(dnormdri,3,3,nmax,"ILPGrapheneHBN:dnormdri");

    memory->create(dnormal,3,3,3,nmax,"ILPGrapheneHBN:dnormal");

  }

  inum = list->inum;

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

   create ILP neighbor list from main neighbor list

   ILP neighbor list stores neighbors of ghost atoms

   ILP_numneigh for calcualting normals and

   ILP_pair_numneigh for calculating force

   create ILP neighbor list from main neighbor list to calcualte normals

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

void PairILPGrapheneHBN::ILP_neigh()

    maxlocal = atom->nmax;

    memory->destroy(ILP_numneigh);

    memory->sfree(ILP_firstneigh);

    memory->create(ILP_numneigh,maxlocal,"ILP:numneigh");

    ILP_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *),"ILP:firstneigh");

    memory->create(ILP_numneigh,maxlocal,"ILPGrapheneHBN:numneigh");

    ILP_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *),"ILPGrapheneHBN:firstneigh");

  }

  allnum = list->inum + list->gnum;

    ILP_firstneigh[i] = neighptr;

    ILP_numneigh[i] = n;

    if (n > 3) error->all(FLERR,"There are too many neighbors for some atoms, please reduce the cutoff for normals");

    ipage->vgot(n);

    if (ipage->status())

      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");