| 查看: 1734 | 回復(fù): 16 | ||
[求助]
最優(yōu)化計(jì)算求助。1stopt或者M(jìn)atlab)
|
|
本人有一個(gè)最優(yōu)化計(jì)算的問題。用1stopt1.5注冊(cè)版編寫的代碼。但是運(yùn)行的時(shí)候總是沒有任何顯示。故意將代碼改錯(cuò),編譯的時(shí)候的確能報(bào)告出錯(cuò)。 另外,貌似1stopt的pascal編程與標(biāo)準(zhǔn)pascal語言有很多不同,比如: 不能自定義結(jié)構(gòu)類型數(shù)據(jù),不能自定義函數(shù),自定義過程貌似也有問題。等等。 請(qǐng)大俠幫忙看看。另外,不知道MatLab能否處理這一問題? 先大致說說這個(gè)問題的物理本質(zhì)和程序的思路。 本質(zhì)就是計(jì)算三個(gè)分子間的范德華力相互作用能。 每個(gè)分子有24個(gè)原子。每個(gè)原子有x,y坐標(biāo)和范德華力相關(guān)參數(shù)兩個(gè)。 所有信息存放在三維數(shù)組里面。數(shù)組的第一指標(biāo)表示不同分子;第二指標(biāo)表示分子內(nèi)各個(gè)原子;第三指標(biāo)表示每個(gè)原子有x,y坐標(biāo)和范德華力相關(guān)參數(shù)兩個(gè)。 M0為分子標(biāo)準(zhǔn)的位置(0,0)和角度。 M1,M2,M3具有相同的分子角度,即將M0的標(biāo)準(zhǔn)角度繞(0,0)點(diǎn)旋轉(zhuǎn)phi角度, M1,M2,M3分別為放置在(0,0), (b,0),和(a*cos(theta),a*sin(theta))位置。 然后分別計(jì)算不同分子各個(gè)原子對(duì)間的范德華力作用能。 程序是用Pascal語言編寫,如下: Title "Type your title here" Parameters a[1,4], b[1,4],theta[0,pi/2],phi[-pi/8,pi/8]; Minimum; StartProgram; const xx: array[1..24] of double =(0.703,1.698,1.698,0.703,-0.703,-1.698,-1.698,-0.703,1.246,3.008,3.008,1.246,-1.246, -3.008,-3.008,-1.246,0,-2.926,-4.138,-2.926,0,2.926,4.138,2.926); yy: array[1..24] of double =(1.698,0.703,-0.703,-1.698,-1.698,-0.703,0.703,1.698,3.008,1.246,-1.246,-3.008,-3.008, -1.246,1.246,3.008,4.138,2.926,0,-2.926,-4.138,-2.926,0,2.926); var i,j: Integer; Etot,r,angle,out1,out2,out3: double; Ax,Ay,Ar,Ae,Bx,By,Br,Be,dist,rou: double; M: array[0..3,1..24,1..4] of double; Begin for i:=1 to 24 do begin M[0,i,1]:=xx; M[0,i,2]:=yy; end; for i:=1 to 8 do begin M[0,i,3]:= 1.9; M[0,i,4]:= 0.044; M[0,i+8,3]:= 1.9; M[0,i+8,4]:= 0.044; M[0,i+16,3]:= 2.11; M[0,i+16,4]:= 0.202; end; M[1]:=M[0]; for i:=1 to 24 do begin r:=sqrt(M[1,i,1]*M[1,i,1]+M[1,i,2]*M[1,i,2]); if M[1,i,1]>=0 then angle:=arcsin(M[1,i,2]/r) else angle:= (M[1,i,2]/abs(M[1,i,2]))*(pi-abs(arcsin(M[1,i,2]/r))); angle:=angle+phi; M[1,i,1]:=r*cos(angle); M[1,i,2]:=r*sin(angle); end; M[2]:=M[1]; M[3]:=M[1]; for i:=1 to 24 do begin M[2,i,1]:=M[1,i,1]+b; M[3,i,1]:=M[1,i,1]+a*cos(theta); M[3,i,2]:=M[1,i,2]+a*sin(theta); end; Etot:=0; for i:=1 to 24 do for j:=1 to 24 do begin Ax:=M[1,i,1]; Ay:=M[1,i,2]; Ar:=M[1,i,3]; Ae:=M[1,i,4]; Bx:=M[2,j,1]; By:=M[2,j,2]; Br:=M[2,j,3]; Be:=M[2,j,4]; dist:=sqrt((Ax-Bx)*(Ax-Bx)+(Ay-By)*(Ay-By)); rou:=dist/(Ar+Br); if dist>3.311 then out1:= sqrt(Ae*Be)*(290000*exp(-12.5*rou)-2.25/exp(6*ln(rou))) else out1:= 336.176*sqrt(Ae*Be)/(rou*rou); Bx:=M[3,j,1]; By:=M[3,j,2]; Br:=M[3,j,3]; Be:=M[3,j,4]; dist:=sqrt((Ax-Bx)*(Ax-Bx)+(Ay-By)*(Ay-By)); rou:=dist/(Ar+Br); if dist>3.311 then out2:= sqrt(Ae*Be)*(290000*exp(-12.5*rou)-2.25/exp(6*ln(rou))) else out2:= 336.176*sqrt(Ae*Be)/(rou*rou); Ax:=M[2,i,1]; Ay:=M[2,i,2]; Ar:=M[2,i,3]; Ae:=M[2,i,4]; dist:=sqrt((Ax-Bx)*(Ax-Bx)+(Ay-By)*(Ay-By)); rou:=dist/(Ar+Br); if dist>3.311 then out3:= sqrt(Ae*Be)*(290000*exp(-12.5*rou)-2.25/exp(6*ln(rou))) else out3:= 336.176*sqrt(Ae*Be)/(rou*rou); Etot:=out1+out2+out3+Etot; end; FunctionResult:= Etot; End; EndProgram; [ Last edited by wdxiao on 2012-3-26 at 14:06 ] |
鐵桿木蟲 (職業(yè)作家)
|
xx,yy數(shù)組分別賦值給M[0,1..24,1]和M[0,1..24,2],表示標(biāo)準(zhǔn)位置的分子內(nèi)各個(gè)原子的x,y坐標(biāo)。對(duì)應(yīng)的表示范德華力的參數(shù)M[0,1..24,3]和M[0,1..24,4]由第一個(gè)for...do循環(huán)賦值。(分子的1~16號(hào)原子是同一種原子,17~24號(hào)原子是另一種原子)。 M0賦值給M1。然后M1旋轉(zhuǎn)phi角度后,賦值給M2和M3。然后M2和M3的各個(gè)原子坐標(biāo)經(jīng)過平移變換。 |
|
xx,yy數(shù)組分別賦值給M[0,1..24,1]和M[0,1..24,2],表示標(biāo)準(zhǔn)位置的分子內(nèi)各個(gè)原子的x,y坐標(biāo)(第一個(gè)for....do循環(huán))。對(duì)應(yīng)的表示范德華力的參數(shù)M[0,1..24,3]和M[0,1..24,4]由第二個(gè)for...do循環(huán)賦值。(分子的1~16號(hào)原子是同一種原子,17~24號(hào)原子是另一種原子)。 M0賦值給M1。然后M1旋轉(zhuǎn)phi角度后,賦值給M2和M3。然后M2和M3的各個(gè)原子坐標(biāo)經(jīng)過平移變換。 |
|
xx,yy數(shù)組分別賦值給M[0,1..24,1]和M[0,1..24,2],表示標(biāo)準(zhǔn)位置的分子內(nèi)各個(gè)原子的x,y坐標(biāo)(第一個(gè)for....do循環(huán))。對(duì)應(yīng)的表示范德華力的參數(shù)M[0,1..24,3]和M[0,1..24,4]由第二個(gè)for...do循環(huán)賦值。(分子的1~16號(hào)原子是同一種原子,17~24號(hào)原子是另一種原子)。 M0賦值給M1。然后M1旋轉(zhuǎn)phi角度后(第三個(gè)for...do循環(huán)),賦值給M2和M3。然后M2和M3的各個(gè)原子坐標(biāo)經(jīng)過平移變換(第四個(gè)for.....do 循環(huán))。 另外,很奇怪的是,我的源程序上第一個(gè)for.....do 循環(huán)為: for i:=1 to 24 do begin M[0,i,1]:=xx; M[0,i,2]:=yy; end; 發(fā)帖的時(shí)候自動(dòng)去掉了,不知道為啥。 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 085410人工智能專碩317求調(diào)劑(0854都可以) +3 | xbxudjdn 2026-03-18 | 3/150 |
|
|---|---|---|---|---|
|
[考研] 一志愿天津大學(xué)化學(xué)工藝專業(yè)(081702)315分求調(diào)劑 +10 | yangfz 2026-03-17 | 10/500 |
|
|
[考研] 收復(fù)試調(diào)劑生 +4 | 雨后秋荷 2026-03-18 | 4/200 |
|
|
[考研] 材料專碩274一志愿陜西師范大學(xué)求調(diào)劑 +6 | 薛云鵬 2026-03-13 | 6/300 |
|
|
[考研] 材料與化工求調(diào)劑 +6 | 為學(xué)666 2026-03-16 | 6/300 |
|
|
[碩博家園] 湖北工業(yè)大學(xué) 生命科學(xué)與健康學(xué)院-課題組招收2026級(jí)食品/生物方向碩士 +3 | 1喜春8 2026-03-17 | 5/250 |
|
|
[考研] 290求調(diào)劑 +3 | p asserby. 2026-03-15 | 4/200 |
|
|
[考研] 材料與化工專碩調(diào)劑 +5 | heming3743 2026-03-16 | 5/250 |
|
|
[考研] 藥學(xué)383 求調(diào)劑 +3 | 藥學(xué)chy 2026-03-15 | 4/200 |
|
|
[考研] 304求調(diào)劑 +5 | 素年祭語 2026-03-15 | 5/250 |
|
|
[考研] 一志愿211 0703方向310分求調(diào)劑 +3 | 努力奮斗112 2026-03-15 | 3/150 |
|
|
[考研] 070300化學(xué)學(xué)碩求調(diào)劑 +6 | 太想進(jìn)步了0608 2026-03-16 | 6/300 |
|
|
[考研] 327求調(diào)劑 +6 | 拾光任染 2026-03-15 | 11/550 |
|
|
[考研] 294求調(diào)劑 +3 | Zys010410@ 2026-03-13 | 4/200 |
|
|
[考研] 288求調(diào)劑 +4 | 奇點(diǎn)0314 2026-03-14 | 4/200 |
|
|
[考研] 080500,材料學(xué)碩302分求調(diào)劑學(xué)校 +4 | 初識(shí)可樂 2026-03-14 | 5/250 |
|
|
[考研] 本科南京大學(xué)一志愿川大藥學(xué)327 +3 | 麥田耕者 2026-03-14 | 3/150 |
|
|
[考研] 復(fù)試調(diào)劑 +4 | z1z2z3879 2026-03-14 | 5/250 |
|
|
[考研] [0860]321分求調(diào)劑,ab區(qū)皆可 +4 | 寶貴熱 2026-03-13 | 4/200 |
|
|
[考研] 材料專碩350 求調(diào)劑 +4 | 王金科 2026-03-12 | 4/200 |
|