| 查看: 1790 | 回复: 16 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[求助]
最优化计算求助！（1stopt或者Matlab）
|
|||
|
本人有一个最优化计算的问题。用1stopt1.5注册版编写的代码。但是运行的时候总是没有任何显示。故意将代码改错，编译的时候的确能报告出错。 另外，貌似1stopt的pascal编程与标准pascal语言有很多不同，比如： 不能自定义结构类型数据，不能自定义函数，自定义过程貌似也有问题。等等。 请大侠帮忙看看。另外，不知道MatLab能否处理这一问题？ 先大致说说这个问题的物理本质和程序的思路。 本质就是计算三个分子间的范德华力相互作用能。 每个分子有24个原子。每个原子有x,y坐标和范德华力相关参数两个。 所有信息存放在三维数组里面。数组的第一指标表示不同分子；第二指标表示分子内各个原子；第三指标表示每个原子有x,y坐标和范德华力相关参数两个。 M0为分子标准的位置（0,0）和角度。 M1,M2,M3具有相同的分子角度，即将M0的标准角度绕(0,0)点旋转phi角度, M1,M2,M3分别为放置在(0,0), (b,0),和(a*cos(theta),a*sin(theta))位置。 然后分别计算不同分子各个原子对间的范德华力作用能。 程序是用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 ] |
回复此楼» 猜你喜欢
重庆交大26年硕士生招生拟调剂通知已出！欢迎加入机器视觉与3D光学成像课题组。
已经有0人回复
-
**
已经有1人回复
-
物理学I论文润色/翻译怎么收费?
已经有248人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有3人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有2人回复
-
B区学生调剂-兰州交通大学材料科学与工程学院
已经有8人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有12人回复
-
桂林理工大学物理学专业招收调剂
已经有11人回复
-
VASP 的一组 GPU / CPU 基准测试记录
已经有0人回复
-
美国顶级物理期刊《应用物理快报》（APL）的编辑欺骗和歧视作者及AIP的官僚主义傲慢
已经有1人回复
» 本主题相关商家推荐: (我也要在这里推广)
|
谢谢！ 这的确是个问题。 我做了相应的修改，如下： 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 if M[1,i,2]>=0 then angle:= pi-abs(arcsin(M[1,i,2]/r)) else angle:= -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); 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; 但运行程序的时候还是没有任何反应。 |
13楼2012-03-26 15:33:57
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
2楼2012-03-26 14:28:57
3楼2012-03-26 14:47:22
4楼2012-03-26 14:50:40
5