| 查看: 1770 | 回复: 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 ] |
回复此楼» 猜你喜欢
070200求调剂，一志愿某211，288分
已经有20人回复
-
法国博士后职位
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有76人回复
-
重庆交大26年硕士生招生拟调剂通知已出！欢迎加入机器视觉与3D光学成像课题组。
已经有0人回复
-
**
已经有1人回复
-
0702一志愿吉大B区求调剂 本科期间发表一篇Sci
已经有3人回复
-
基底STO,薄膜SRO，XRD里面的振荡，是laue震荡还是kiessig振荡? 怎么判断？
已经有2人回复
-
B区学生调剂-兰州交通大学材料科学与工程学院
已经有8人回复
-
山西大同大学物理学专业还有调剂名额，欢迎调剂！
已经有9人回复
-
桂林理工大学物理学专业招收调剂
已经有6人回复
» 本主题相关商家推荐: (我也要在这里推广)
» 本主题相关价值贴推荐，对您同样有帮助:
- 1stOpt求解常微分方程边值问题
已经有5人回复
- 哪位大哥能帮我用1stopt3.0以上版本跑一个拟合程序
已经有11人回复
- 1stOpt的函数作图不靠谱吗？
已经有4人回复
- 求助有关Matlab有约束非线性最优化问题
已经有9人回复
- 用origin或1stopt拟合一条隐函数方程曲线
已经有4人回复
- 【求助】用matlab最优化方法进行参数拟合
已经有17人回复
- 求助关于1stopt的问题
已经有6人回复
- 【求助】MUSIC的MD模块计算甲烷在IRMOF-1上的扩散系数
已经有11人回复
- 用1stopt进行曲线拟合 每次给出结果不一样
已经有6人回复
- 用1stopt进行数据拟合时，无法运行sharedmodel ,求助高手
已经有7人回复
- 求助关于1stopt拟合~··~~
已经有4人回复
- 【求助】精通matlab最优化计算，一书中的funval函数是什么意思
已经有10人回复
- 【求助】如何计算有机小分子的三重态能级（T1-S0）？
已经有12人回复
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
|
xx,yy数组分别赋值给M[0,1..24,1]和M[0,1..24,2]，表示标准位置的分子内各个原子的x,y坐标（第一个for....do循环）。对应的表示范德华力的参数M[0,1..24,3]和M[0,1..24,4]由第二个for...do循环赋值。（分子的1~16号原子是同一种原子，17~24号原子是另一种原子）。 M0赋值给M1。然后M1旋转phi角度后（第三个for...do循环），赋值给M2和M3。然后M2和M3的各个原子坐标经过平移变换（第四个for.....do 循环）。 另外，很奇怪的是，我的源程序上第一个for.....do 循环为： for i:=1 to 24 do begin M[0,i,1]:=xx; M[0,i,2]:=yy; end; 发帖的时候自动去掉了，不知道为啥。 |
5楼2012-03-26 14:56:46
6楼2012-03-26 15:10:32
7楼2012-03-26 15:11:09
8楼2012-03-26 15:11:36
9楼2012-03-26 15:12:17
10楼2012-03-26 15:12:53
5