24小时热门版块排行榜    

查看: 1810  |  回复: 7

frcon

新虫 (初入文坛)

[求助] LAMMPS源代码的求助 已有1人参与

给位老师,大家好,最近我在修改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]的意思,各位有知道的能讲解一下吗?万分感谢!
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lsloneil

专家顾问 (正式写手)

【答案】应助回帖

感谢参与,应助指数 +1
这很好理解啊,模拟中可以定义多种类型的键,每一种类型的键都有相应的力常数k和平衡距离r0, type就表示健的种类。

p.s. 如果不懂c++,还是补一点基础知识再去看源代码。
2楼2016-02-23 15:27:18
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

frcon

新虫 (初入文坛)

引用回帖:
2楼: Originally posted by lsloneil at 2016-02-23 15:27:18
这很好理解啊,模拟中可以定义多种类型的键,每一种类型的键都有相应的力常数k和平衡距离r0, type就表示健的种类。

p.s. 如果不懂c++,还是补一点基础知识再去看源代码。

好的,谢谢
3楼2016-02-23 15:33:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

frcon

新虫 (初入文坛)

引用回帖:
2楼: Originally posted by lsloneil at 2016-02-23 15:27:18
这很好理解啊,模拟中可以定义多种类型的键,每一种类型的键都有相应的力常数k和平衡距离r0, type就表示健的种类。

p.s. 如果不懂c++,还是补一点基础知识再去看源代码。

那像代码里的ebond,fbond这些定义对应的意思可以在哪里找到
4楼2016-02-23 15:35:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lsloneil

专家顾问 (正式写手)

【答案】应助回帖

引用回帖:
4楼: Originally posted by frcon at 2016-02-22 19:35:20
那像代码里的ebond,fbond这些定义对应的意思可以在哪里找到...

没法找到,lammps代码不可能每一行都给你注释,每个变量的具体含义只能根据变量的名字和在程序中的调用,结合注释才能推测。

像ebond, fbond都还是比较简单的能根据名字能推测出含义的变量,也就是键的能量和力。如果碰到其他你看了半天代码也推测不出来含义的变量,可以到lammps mailing list去请教别人。lammps developer guide(http://lammps.sandia.gov/doc/Developer.pdf)可以帮助你了解代码的结构,但不会告诉你每个变量每个函数的具体含义。总之要了解代码是要花一些时间的。
5楼2016-02-24 00:26:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

frcon

新虫 (初入文坛)

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

发自小木虫IOS客户端
6楼2016-02-24 00:30:42
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lsloneil

专家顾问 (正式写手)

【答案】应助回帖

引用回帖:
6楼: Originally posted by frcon at 2016-02-23 04:30:42
好的,谢谢你,我已经能把代码和势函数对应上了,还有个疑问麻烦了,就是力的公司是怎么写进去的,因为我不是学习分子动力学的,所以不是很清楚,了解的还不够

这很简单,只要你知道了你键能还有受力(也就是键能对距离的导数取负值),直接把ebond和fbond的函数形式改掉就行了。以bond_harmonic.cpp为例,在compute()这个函数里,你可以看到以下代码
CODE:
    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之间的区别,这个不难的。
7楼2016-02-24 03:22:14
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

frcon

新虫 (初入文坛)

引用回帖:
7楼: Originally posted by lsloneil at 2016-02-24 03:22:14
这很简单,只要你知道了你键能还有受力(也就是键能对距离的导数取负值),直接把ebond和fbond的函数形式改掉就行了。以bond_harmonic.cpp为例,在compute()这个函数里,你可以看到以下代码
    rsq = delx*delx  ...

谢谢你,我比较了这两个代码,发现了他们的不同,也知道了怎么去修改,非常感谢,我去试试

发自小木虫IOS客户端
8楼2016-02-24 15:16:59
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 frcon 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 Vi50GxzrFcSG 2026-08-07 6/300 2026-08-08 10:19 by 3OOjAIS77qg2
[基金申请] fileCode有新解读? +8 Tide man 2026-08-08 16/800 2026-08-08 09:59 by Tide man
[硕博家园] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 HEQlVqMTIA7d 2026-08-07 5/250 2026-08-08 08:59 by 6vVgjDL4CnGu
[基金申请] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 KXLV3nuBVBY7 2026-08-07 4/200 2026-08-08 08:21 by 6vVgjDL4CnGu
[基金申请] 基金中了 +14 laoda193707 2026-08-06 14/700 2026-08-08 00:23 by 实验小白ha
[基金申请] 关于filecode +4 布布和一二 2026-08-07 7/350 2026-08-07 22:55 by zhanghaozhu
[基金申请] 娱乐 +6 Tide man 2026-08-03 6/300 2026-08-07 22:40 by 铁帽子农民
[基金申请] 化学口download_prp&fileCode的固定段好像这几天一直没变,有变的大神么? +3 Tide man 2026-08-07 4/200 2026-08-07 22:39 by Tide man
[基金申请] 固定端突然变了,今天 +6 archvillain 2026-08-06 10/500 2026-08-07 16:03 by 医学老男孩
[基金申请] filecode变化情况 +6 布布和一二 2026-08-07 22/1100 2026-08-07 14:45 by 且听虎啸
[论文投稿] 十年后又回来了,论文投稿求助 +4 哈哈114477 2026-08-01 4/200 2026-08-07 14:39 by jgy194592
[基金申请] 大家散了吧,后缀研究没有意义,别浪费时间了,过好目前的每一天,不要焦虑 +5 Tide man 2026-08-06 7/350 2026-08-07 13:11 by 医学老男孩
[基金申请] filecode +14 等待解的谜 2026-08-06 19/950 2026-08-07 12:20 by wlwhappy
[基金申请] filecode +8 布布和一二 2026-08-06 11/550 2026-08-06 20:41 by tangpu318
[基金申请] 求各位大神看下 100+6 hpkpkpkp 2026-08-05 33/1650 2026-08-06 14:49 by zhiyanjiang
[基金申请] 影响面上的因素 +8 布布和一二 2026-08-05 11/550 2026-08-06 10:41 by 宝贝虫子
[基金申请] 有没有H口的?有收到消息的吗? +3 超级海虾 2026-08-04 3/150 2026-08-04 17:26 by 学教育滴
[基金申请] 纯娱乐,不喜欢勿喷 +7 Tide man 2026-08-04 10/500 2026-08-04 15:10 by loufangrui
[基金申请] 面上提前没消息,有中的吗 +14 archvillain 2026-08-02 18/900 2026-08-04 14:42 by archvillain
[基金申请] 什么时候能放榜呀? +3 Jacob678 2026-08-03 3/150 2026-08-03 16:14 by gltch
信息提示
请填处理意见