| 查看: 2790 | 回復: 15 | |||
[交流]
【求助】使用MATLAB怎么實現(xiàn)擬合力場參數(shù)的程序?【已解決】 已有6人參與
|
|
可以將根據(jù)VDW公式計算出的E_VDW_total視為各個原子間距離r(r是向量)以及原子VDW參數(shù){sigma}、{epsilon}變量集的函數(shù)。r (i)代表第i個掃描點的原子間距離向量,定義誤差函數(shù)ErrF=∑(r(i)下量化計算的能量-r(i)下的E_VDW_total)^2,然后對ErrF對{sigma}和{epsilon}變量集中的每個變量求導得0,如果寫出來是線性形式可以用矩陣方程來解,如果不是的話可以用非線性優(yōu)化方法解,解出{sigma}、{epsilon},用這樣參數(shù)計算的VDW作用可以對所有掃描的點的能量都能較好描述。可能說得比較抽象,這與擬合ESP電荷的方法比較像,可看相關文獻。 [ Last edited by nono2009 on 2010-12-1 at 08:20 ] |


木蟲 (小有名氣)
Matlab

|
我是在問問題,這個是我寫的程序,其實跟最小二乘法有關: function [a1,b1,a2,b2]=forcefieldfitting(r,e) %使用DFT計算的數(shù)據(jù)擬合力場參數(shù),e是使用量化計算出現(xiàn)的相互作用能,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的維數(shù)不相等!'); return; end %維數(shù)檢查 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); |

| 但是就是運行不出來,我是先用量子化學計算相互作用能,我要計算的是CO2的LJ參數(shù)勢,計算出來的能量使用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)求導,然后命令這些導數(shù)等于0,然后利用矩陣的左除,不知道這么做算法對不對? |



金蟲 (職業(yè)作家)



| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 290分調劑求助 +6 | 吉祥止止陳 2026-03-25 | 6/300 |
|
|---|---|---|---|---|
|
[考研] 070300化學354求調劑 +12 | 101次希望 2026-03-28 | 12/600 |
|
|
[考研] 08工科,295,接受跨專業(yè)調劑 +3 | lmnlzy 2026-03-31 | 3/150 |
|
|
[考研] 311求調劑一志愿合肥工業(yè)大學 +8 | 秋二十二 2026-03-30 | 8/400 |
|
|
[考研] 085602 307分 求調劑 +8 | 不知道叫什么! 2026-03-26 | 8/400 |
|
|
[論文投稿]
chinese chemical letters英文版投稿求助
130+3
|
Yishengeryi 2026-03-30 | 3/150 |
|
|
[考研] 08工科,295,接受跨專業(yè)調劑 +6 | lmnlzy 2026-03-30 | 6/300 |
|
|
[考研] 085601一志愿西北工業(yè)大學初試346 +4 | 085601初試346 2026-03-30 | 4/200 |
|
|
[考研] 哈爾濱工業(yè)大學材料與化工專碩378求調劑 +3 | 塔比烏斯 2026-03-30 | 3/150 |
|
|
[考研] 抱歉 +4 | 田洪有 2026-03-30 | 4/200 |
|
|
[考研] 328求調劑 +8 | 嗯滴的基本都 2026-03-27 | 8/400 |
|
|
[考研] 322求調劑:一志愿湖南大學 材料與化工(085600),已過六級。 +9 | XX小鄧 2026-03-29 | 9/450 |
|
|
[考研] 085602 化學工程專碩 340分求調劑 +4 | qianbai11 2026-03-29 | 4/200 |
|
|
[考研] 材料與化工272求調劑 +21 | 阿斯蒂芬2004 2026-03-28 | 21/1050 |
|
|
[考研] 調劑考研 +3 | 王杰一 2026-03-29 | 3/150 |
|
|
[考研] 081200-11408-276學碩求調劑 +6 | 崔wj 2026-03-26 | 6/300 |
|
|
[考研] 299求調劑 +7 | 嗯嗯嗯嗯2 2026-03-27 | 7/350 |
|
|
[考研] 一志愿南京航空航天大學材料學碩求調劑 +3 | @taotao 2026-03-28 | 3/150 |
|
|
[考研] 一志愿南師大0703化學 275求調劑 +4 | Ripcord上岸 2026-03-27 | 4/200 |
|
|
[考研] 286求調劑 +4 | lim0922 2026-03-26 | 4/200 |
|