| 查看: 2803 | 回复: 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 ] |
回复此楼» 猜你喜欢
复试调剂
已经有7人回复
-
一志愿北京化工085600 310分求调剂
已经有13人回复
-
085600，专业课化工原理，320分求调剂
已经有10人回复
-
288求调剂
已经有5人回复
-
0703化学321分求调剂
已经有16人回复
-
318求调剂
已经有13人回复
-
一志愿南航，数一英一学硕317求调剂！！
已经有5人回复
-
一志愿9材料学硕297已过六级求调剂推荐
已经有11人回复
-
085602调剂 初试总分335
已经有5人回复
-
302分 085601求调剂推荐
已经有8人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- matlab拟合拟均相动力学参数
已经有6人回复
- 需要用matlab拟合数据的外行人急需帮助
已经有10人回复
- 求助啊！求一程序，用matlab程序做，用蒙特卡洛方法模拟
已经有10人回复
- matlab非线性参数拟合问题
已经有7人回复
- 如何用MATLAB 实现化学反应方程式（写程序代码）?
已经有20人回复
- Matlab人工神经网络工具箱的位置及使用
已经有5人回复
- 看不懂matlab程序该怎么办？
已经有7人回复
- 【求助】如何根据势能曲线拟合力场参数
已经有19人回复
2楼2010-05-31 19:05:46
lijinfeng042
木虫 (小有名气)
Matlab
- 仿真EPI: 2
- 应助: 1 (幼儿园)
- 金币: 2156.1
- 散金: 115
- 帖子: 291
- 在线: 31.5小时
- 虫号: 1019062
- 注册: 2010-05-15
- 性别: GG
- 专业: 通信理论与系统
3楼2010-05-31 20:01:07
|
我是在问问题，这个是我写的程序，其实跟最小二乘法有关： 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
6楼2010-06-01 15:00:55
7楼2010-06-01 15:01:36
onesupeng
金虫 (职业作家)
- 仿真EPI: 11
- 应助: 256 (大学生)
- 贵宾: 1.36
- 金币: 2336.2
- 散金: 9228
- 红花: 92
- 帖子: 4583
- 在线: 1303.8小时
- 虫号: 394701
- 注册: 2007-06-07
- 专业: 流体力学
8楼2010-06-02 06:56:19
9楼2010-06-02 08:19:23
10楼2010-06-02 16:36:13
