24小时热门版块排行榜    

查看: 1857  |  回复: 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一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 ASdOkHsho7FD 2026-08-28 6/300 2026-08-29 07:59 by gy1nBQXYQJqL
[基金申请] 系统查不到 +11 董八千 2026-08-26 11/550 2026-08-28 18:06 by Leogzhya
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +8 苏东坡二世 2026-08-23 8/400 2026-08-28 14:41 by lkforward
[基金申请] 基金不中,共勉 +11 eulota 2026-08-26 11/550 2026-08-28 14:22 by 火星超人xi
[基金申请] 麻烦专家们看看评委们的意见(F口面上) +4 gdd2018 2026-08-28 9/450 2026-08-28 14:22 by jklily
[基金申请] 哪位高人中了,把查询到的截图贴出来让我看看,让我长长见识 +5 yuleib84 2026-08-26 6/300 2026-08-28 00:02 by yudaoqian88
[基金申请] 面上合作单位盖章 +5 ssyjh 2026-08-27 5/250 2026-08-27 20:50 by gdfollow
[基金申请] 申请删除本帖 +6 lyz123lyz 2026-08-27 7/350 2026-08-27 17:31 by 宁静致远sy
[基金申请] 基金未中,这种答复是模板吗? +5 zhaosm1982 2026-08-27 6/300 2026-08-27 16:00 by lfy8008
[文学芳草园] 梦想 +3 myrtle 2026-08-26 3/150 2026-08-27 10:01 by angelyueyi
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件,有关系吗? +25 木水思豆 2026-08-25 28/1400 2026-08-26 14:53 by draco1987
[基金申请] 为什么国自然不能直接公布 +4 bjdxyxy 2026-08-26 4/200 2026-08-26 13:12 by qingmu1201
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 国际合作可查了,中了面上 (EPI+1)(金币+50) +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 国合现在查不到了吗? +10 chengyan1220 2026-08-24 20/1000 2026-08-26 08:57 by peasantsprig
[基金申请] 在坚冰还盖着北海的时候,我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 如果此刻你正在为国基感到焦虑,不妨来听听这首《基金之外》 +8 scalable 2026-08-24 8/400 2026-08-25 12:52 by jnhyjjm
信息提示
请填处理意见