LAMMPS源代码- 源自精小木虫论坛华贴

Posted 2020-10-30 Simyang

关于源代码

Ask：想看看lammps源代码里面是如何计算原子间相互作用的，结果看来看去都没发现哪块是计算相互作用的，lammps的源代码写的很不容易看清楚啊。 从主程序main.cpp出发，里面除去MPI相关的函数，和创建一个lammps实例，剩下就一句话  lammps->input->file();   调用lammps下的input类的file函数，处理输入文件。 这个file函数里面除去很多判断最后就落到parse函数，对命令进行解析，然后进入execute_command函数调用相应的命令处理函数，但是这里面没有run命令的处理函数。 以pair_style命令为例，对应的pair_style函数会将输入的命令参数（力场名称）和类force下的pair_style字符串比较（然而这个字符串的赋值始终没找到），如果符合就进入相应的settings函数。然后，这个settings函数里面一般就一行判断参数个数的语句就没了。 求高手指点一下，lammps源代码的结构是怎么样的？哪个文件里面包含了原子间相互作用的计算？

Response1：用的是C++的多态技术；

Response2：比如，想看lj96势函数的，找pair_lj98_cut.cpp，打开这个CPP文件，就能看到计算作用力的代码；

Response3：嗯，这个应该是可以理解的，不过我想看的是lammps中怎么把计算的力结果返回到主程序的，分子动力学有个力->速度->位置->力的循环迭代，应该是在主程序中完成的吧，我想知道这个用不同势函数计算的力结果和主程序之间的接口。

Response4：也许是在积分的函数里面调用的。积分有两步，中间夹有求力。我自己写程序时就是将求力的函数交给积分的函数调用的。当然，你如果只看主函数可能看不明白，这种大型程序一般会包装好几层才会到具体的真正做计算的函数。我还没仔细查看lammps代 码，以上只是我的猜测。

Response5：楼主应该是看过一些源代码了吧，程序从main函数开始，然后是就转到类 input, 在这个类中，完成了文件的读取，并解析输入命令，然后这个类中的成员函数， excute_command()就是执行每个模拟指令的，计算力的过程也是在这里完成。input类继承了   Pointers类，其中的类指钍Force *force可以调用不同的pair_style, bond_style类，当然对应的计算作用力的程序也就在这个过程中完成了  看懂源代码需要花一翻功夫的，我也只看粗略看了一下，希望对你有帮助

Response6：嗯，今天又看了一些代码，发现了run函数，可能是具体的计算过程由update类下面的integrate类的run函数完成，但是这个run函数我目前看到的是虚函数，还没发现可以实例化的代码。然而关于主程序还是不太明白，run函数是在实例化lammps类的过程中完成调用的，但是之后才有input类下的file函数关于输入文件的处理，但是，file函数之后就直接delete lammps，删除实例了！这样的话真正的运算在哪里？

LAMMPS源代码的求助

Ask：

最近我在修改lammps中bond_harmonic的源代码，对其中部分不是很清楚，在这里请教一下：

bond_harmonic的势函数为E=K（r-r0)^2, 在源代码中

rsq = delx*delx + dely*dely + delz*delz;

r = sqrt(rsq);

dr = r - r0[type]; rk = k[type] * dr;

其势函数是这样写的，在这里，我不是很懂程序中[type]的意思，各位有知道的能讲解一下吗？万分感谢！

Response1：这很好理解啊，模拟中可以定义多种类型的键，每一种类型的键都有相应的力常数k和平衡距离r0， type就表示健的种类。

Response2：

那像代码里的ebond，fbond这些定义对应的意思可以在哪里找到...

没法找到，lammps代码不可能每一行都给你注释，每个变量的具体含义只能根据变量的名字和在程序中的调用，结合注释才能推测。像ebond， fbond都还是比较简单的能根据名字能推测出含义的变量，也就是键的能量和力。如果碰到其他你看了半天代码也推测不出来含义的变量，可以到lammps mailing list去请教别人。lammps developer guide（http://lammps.sandia.gov/doc/Developer.pdf）可以帮助你了解代码的结构，但不会告诉你每个变量每个函数的具体含义。总之要了解代码是要花一些时间的。

Response3：

好的，谢谢你，我已经能把代码和势函数对应上了，还有个疑问麻烦了，就是力的公式是怎么写进去的，因为我不是学习分子动力学的，所以不是很清楚，了解的还不够

 这很简单，只要你知道了你键能还有受力（也就是键能对距离的导数取负值），直接把ebond和fbond的函数形式改掉就行了。以bond_harmonic.cpp为例，在compute()这个函数里，你可以看到以下代码

    rsq = delx*delx + dely*dely + delz*delz;
    r = sqrt(rsq);
    dr = r - r0[type];
    rk = k[type] * dr;

    // force & energy

    if (r > 0.0) fbond = -2.0*rk/r;
    else fbond = 0.0;

    if (eflag) ebond = rk*dr;

 你需要的就是在compute()还有single()两个函数里把fbond和ebond的函数形式改成你所需要的形式。如果你不太确定怎么改，可以参考bond_morse.cpp，也就是把harmonic bond换成morse bond，看看代码是怎么变动的。
 另外如果你定义键能的参数也发生了变动（比方说你需要更多的参数来定义键能），同时也需要对coeff()这个函数进行修改来改变参数传递，具体例子参考bond_morse.cpp和bond_harmonic.cpp之间的区别，这个不难的