比如说我要加一个新的势能函数，那些如何计算力的方法，计算位移的方法我都要重写吗？还是说我只要改一部分就行了。
问题是我不知道内部是如何算得啊，大家有做过这个的吗？

当然不用！

percysong同学我不知道你有没有用过ABAQUS的用户子程序(UMAT,UEL)，是不是就像ABAQUS那样只把几个指定的需要更新的变量在子程序里计算一下，然后返回给主程序就可以了，至于用什么样的算法求解方程组的根我们根本不需要知道，LAMMPS是这样工作的吗？

in case of EAM
The EAM potentials LAMMPS uses are tabulated in files,
so in principal you can use any functional form for any
material, tabulate it, and read it in from a file.  The
pair_style eam doc page lists some WWW sites that
have EAM files for different materials in the right
format, or you could convert them.

It you want an analytic EAM form (e.g. Sutton/Chen),
you could implement it, e.g. as pair_style eam/suttonchen
and hard code the formulas.  This is not hard to do;
most pair potentials in LAMMPS are in that mode.

Start with a doc page thta defines the needed inputs
and code up the formulas and send it to me.  It would
then be easy to make a new pair style file.

You'll need to write a doc page similar to doc/pair_eam.txt,
call it pair_eam_spline.txt.  This will force you to define
exactly what parameters you want to input, so it's
helpful in defining what needs to be coded.

potential itself, yes code the energy and force in C
or C++ for the  equations.

steve

lammps只能运用表列好了的势能函数计算，而势能函数的参数，是第一性原理计算拟合出来的，你先要有第一性计算的拟合结果然后按要求写好适合lammps的文件，如果势能函数有解析形式，那就更方便一些，直接编程自己算结果然后加上去就好了。lammps只提供高效的MD算法框架，并不计算生成势能函数。

比如说我在SRC目录下发现了pair_tersoff.cpp文件和pair_tersoff.h文件，如果我想用一个是能函数比如说叫做bang势能函数，那我就要写两个文件一个是pair_bang.cpp另一个是pair_bang.h对吗？
其中对照pair_tersoff.cpp来改写pair_bang.cpp， 对照pair_tersoff.h来改写pair_bang.h对吗？
只要写完了这两个文件就可以了对吗？其他的任何东西都不要动对吗？

[ 本帖最后由 sunysb 于 2008-8-4 10:09 编辑 ]

YES, although you do not change algorithm and structure in force. h, comm.h, memory.h etc. ,you still need to read them though to ensure  the potential you generate correct.

percysong能否写点修改code的经验啊，很希望你能写点

记得我以前问过sutton chen的问题，目前结果不好，skin部分的位势截断和源代码不太一样，等有了好结果希望能写点对大家有用的东西。

我看了类似的程序，就是不知道要更新哪些变量，要保留哪些变量。

你具体要加哪个势函数？

我要加Tersoff-Brenner，在SRC下面已经有了Tersoff但是我想要Tersoff-Brenner

我自己正在试着把一个eam势加上去，希望能完善到可以给steve发过去

void PairLJCut::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,r2inv,r6inv,forcelj,factor_lj;
  int *ilist,*jlist,*numneigh,**firstneigh;

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

  double **x = atom->x;
  double **f = atom->f;
  int *type = atom->type;
  int nlocal = atom->nlocal;
  int nall = nlocal + atom->nghost;
  double *special_lj = force->special_lj;
  int newton_pair = force->newton_pair;

  inum = list->inum;
  ilist = list->ilist;
  numneigh = list->numneigh;
  firstneigh = list->firstneigh;
  
  // loop over neighbors of my atoms

  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    xtmp = x[0];
    ytmp = x[1];
    ztmp = x[2];
    itype = type;
    jlist = firstneigh;
    jnum = numneigh;

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

      if (j < nall) factor_lj = 1.0;
      else {
        factor_lj = special_lj[j/nall];
        j %= nall;
      }

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

      if (rsq < cutsq[itype][jtype]) {
        r2inv = 1.0/rsq;
        r6inv = r2inv*r2inv*r2inv;
        forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
        fpair = factor_lj*forcelj*r2inv;

        f[0] += delx*fpair;
        f[1] += dely*fpair;
        f[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 = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
            offset[itype][jtype];
          evdwl *= factor_lj;
        }

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

  if (vflag_fdotr) virial_compute();
}
以上的程序出自src/下的 Pair_lj_cut.cpp, 其中第9行的 if (eflag || vflag) ev_setup(eflag,vflag); 当中eflag和vflag还有ev_setup是什么意思啊？

pair.cpp and pair.h 中可以找到答案

/* ----------------------------------------------------------------------
   setup for energy, virial computation
   see integrate::ev_set() for values of eflag (0-3) and vflag (0-6)
------------------------------------------------------------------------- */

  int evflag;                          // energy,virial settings