| 查看: 3354 | 回复: 15 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[交流]
【求助】使用MATLAB怎么实现拟合力场参数的程序？【已解决】 已有6人参与
|
|||
|
可以将根据VDW公式计算出的E_VDW_total视为各个原子间距离r（r是向量）以及原子VDW参数{sigma}、{epsilon}变量集的函数。r (i)代表第i个扫描点的原子间距离向量，定义误差函数ErrF=∑(r(i)下量化计算的能量-r(i)下的E_VDW_total)^2，然后对ErrF对{sigma}和{epsilon}变量集中的每个变量求导得0，如果写出来是线性形式可以用矩阵方程来解，如果不是的话可以用非线性优化方法解，解出{sigma}、{epsilon}，用这样参数计算的VDW作用可以对所有扫描的点的能量都能较好描述。可能说得比较抽象，这与拟合ESP电荷的方法比较像，可看相关文献。 [ Last edited by nono2009 on 2010-12-1 at 08:20 ] |
回复此楼» 猜你喜欢
今天务委会开完了，明天出结果吗
已经有22人回复
-
售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急
已经有3人回复
-
2026年8月25日国自然放榜前突然收到列入评审专家邮件，有关系吗？
已经有19人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.5.4,科目全,可+急
已经有3人回复
-
在坚冰还盖着北海的时候，我看到了怒放的梅花。
已经有9人回复
-
明天应该可查了！？
已经有6人回复
-
有没有大神帮我看看基金代码
已经有10人回复
-
放榜前的不淡定
已经有19人回复
-
某些机构，以效率低为荣，以效率低作为存在感
已经有10人回复
-
2026国自然函评费到账
已经有23人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- matlab拟合拟均相动力学参数
已经有6人回复
- 需要用matlab拟合数据的外行人急需帮助
已经有10人回复
- 求助啊！求一程序，用matlab程序做，用蒙特卡洛方法模拟
已经有10人回复
- matlab非线性参数拟合问题
已经有7人回复
- 如何用MATLAB 实现化学反应方程式（写程序代码）?
已经有20人回复
- Matlab人工神经网络工具箱的位置及使用
已经有5人回复
- 看不懂matlab程序该怎么办？
已经有7人回复
- 【求助】如何根据势能曲线拟合力场参数
已经有19人回复
lijinfeng042
木虫 (小有名气)
Matlab
- 仿真EPI: 2
- 应助: 1 (幼儿园)
- 金币: 2156.1
- 散金: 115
- 帖子: 291
- 在线: 31.5小时
- 虫号: 1019062
- 注册: 2010-05-15
- 性别: GG
- 专业: 通信理论与系统
3楼2010-05-31 20:01:07
2楼2010-05-31 19:05:46
|
我是在问问题，这个是我写的程序，其实跟最小二乘法有关： function [a1,b1,a2,b2]=forcefieldfitting(r,e) %使用DFT计算的数据拟合力场参数,e是使用量化计算出现的相互作用能，r是两个分子片段之间的距离 %LJ-12-06公式中的EPSILON(C):a1 %LJ-12-06公式中的SIGMA(C):b1 %LJ-12-06公式中的EPSILON(O):a2 %LJ-12-06公式中的SIGMA(O):b2 %单个的CO2的能量可以认为是两个O的能量和一个C的能量的加和 if(length(e)==length(r)) n=length(e); else disp('e和r的维数不相等！'); return; end %维数检查 A=zeros(4,3) B=zeros(4,2); for i=1:n A(1,1)=A(1,1)+(384*a1^2*b1^23)/r(i)^24; A(1,2)=A(1,2)+(192*a1^2*b1^11)/r(i)^12; A(1,3)=A(1,3)+(48*a1*e*b1^5)/r(i)^6; A(2,1)=A(2,1)+(32*a1*b1^24)/r(i)^24; A(2,2)=A(2,2)+(32*a1*b1^12)/r(i)^12 A(2,3)=A(2,3)+(8*e(i)*b1^6)/r(i)^6; A(3,1)=A(3,1)+(1536*a2^2*b2^23)/r(i)^24; A(3,2)=A(3,2)+(768*a2^2*b2^23)/r(i)^12; A(3,3)=A(3,3)+(96*a2*e(i)*b2^5)/r(i)^6; A(4,1)=A(4,1)+(128*a2*b2^24)/r(i)^24; A(4,2)=A(4,2)+(128*a2*b2^12)/r(i)^12; A(4,3)=A(4,3)+(16*e(i)*b2^6)/r(i)^6; B(1,1)=B(1,1)+(576*a1^2*b1^17)/r(i)^18; B(1,2)=B(1,2)+(96*a1*e(i)*b1^11)/r(i)^12; B(2,1)=B(2,1)+(64*a1*b1^18)/r(i)^18; B(2,2)=B(2,2)+(8*e(i)*a1^12)/r(i)^12; B(3,1)=B(3,1)+(2304*a2^2*b2^17)/r(i)^18; B(3,2)=B(3,2)+(192*a2*e(i)*b2^11)/r(i)^12; B(4,1)=B(4,1)+(256*a2*b2^18)/r(i)^18; B(4,3)=B(4,3)+(16*e(i)*b2^12)/r(i)^12; end s=A\B; a1=s(1); b1=s(2); a2=s(3); b2=s(4); |
4楼2010-05-31 22:50:07
★
zzuwangshilei(金币+1):鼓励深入讨论自己的帖子 2010-06-01 14:21:17
zzuwangshilei(金币+1):鼓励深入讨论自己的帖子 2010-06-01 14:21:17
| 但是就是运行不出来，我是先用量子化学计算相互作用能，我要计算的是CO2的LJ参数势，计算出来的能量使用E=4*EPSILON(C)((SIGMA(C)**12/R**12-(SIGMA(C)**6/R**6)+2*EPSILON(O)((SIGMA(O)**12/R**12-(SIGMA(O)**6/R**6)然后把此式分别对EPSILON(C),SIGMA(C),EPSILON(O),SIGMA(O)求导，然后命令这些导数等于0，然后利用矩阵的左除，不知道这么做算法对不对？ |
5楼2010-05-31 22:53:50
