24小时热门版块排行榜    

查看: 1737  |  回复: 16
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

wdxiao

银虫 (正式写手)

[求助] 最优化计算求助!(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 ]
回复此楼

» 猜你喜欢

» 本主题相关商家推荐: (我也要在这里推广)

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wdxiao

银虫 (正式写手)

引用回帖:
2楼: Originally posted by dingd at 2012-03-26 14:28:57:
不知道你的M值是在哪赋值的?

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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 17 个回答

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +1
不知道你的M值是在哪赋值的?
2楼2012-03-26 14:28:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wdxiao

银虫 (正式写手)

引用回帖:
2楼: Originally posted by dingd at 2012-03-26 14:28:57:
不知道你的M值是在哪赋值的?

xx,yy数组分别赋值给M[0,1..24,1]和M[0,1..24,2],表示标准位置的分子内各个原子的x,y坐标。对应的表示范德华力的参数M[0,1..24,3]和M[0,1..24,4]由第一个for...do循环赋值。(分子的1~16号原子是同一种原子,17~24号原子是另一种原子)。

M0赋值给M1。然后M1旋转phi角度后,赋值给M2和M3。然后M2和M3的各个原子坐标经过平移变换。
3楼2012-03-26 14:47:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wdxiao

银虫 (正式写手)

引用回帖:
2楼: Originally posted by dingd at 2012-03-26 14:28:57:
不知道你的M值是在哪赋值的?

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角度后,赋值给M2和M3。然后M2和M3的各个原子坐标经过平移变换。
4楼2012-03-26 14:50:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 求材料调剂 +8 隔壁陈先生 2026-03-12 8/400 2026-03-18 22:19 by li123456789.
[考研] 321求调剂 +3 何润采123 2026-03-18 3/150 2026-03-18 21:27 by li123456789.
[考研] 材料专硕英一数二306 +4 z1z2z3879 2026-03-18 4/200 2026-03-18 20:00 by 楤哥
[考研] 26调剂/材料/英一数二/总分289/已过A区线 +7 步川酷紫123 2026-03-13 7/350 2026-03-18 17:12 by 尽舜尧1
[考研] 085600材料与化工 +5 安全上岸! 2026-03-16 5/250 2026-03-18 15:33 by cmz0325
[考研] 化学工程321分求调剂 +15 大米饭! 2026-03-15 18/900 2026-03-18 14:52 by haxia
[考研] 311求调剂 +6 26研0 2026-03-15 6/300 2026-03-18 14:43 by haxia
[考研] 0817 化学工程 299分求调剂 有科研经历 有二区文章 +7 rare12345 2026-03-18 7/350 2026-03-18 14:31 by laoshidan
[考研] 0854,计算机类招收调剂 +3 胡辣汤放糖 2026-03-15 6/300 2026-03-18 12:09 by 上岸上岸……..
[考研] 0703化学求调剂 总分331 +3 ZY-05 2026-03-13 3/150 2026-03-18 10:58 by macy2011
[考研] 26考研求调剂 +6 丶宏Sir 2026-03-13 6/300 2026-03-17 16:13 by 醉在风里
[论文投稿] 有没有大佬发小论文能带我个二作 +3 增锐漏人 2026-03-17 4/200 2026-03-17 09:26 by xs74101122
[考研] 070303 总分349求调剂 +3 LJY9966 2026-03-15 5/250 2026-03-16 14:24 by xwxstudy
[考研] 277材料科学与工程080500求调剂 +3 自由煎饼果子 2026-03-16 3/150 2026-03-16 14:10 by 运气yunqi
[考研] 22408总分284求调剂 +3 InAspic 2026-03-13 3/150 2026-03-15 11:10 by zhq0425
[考研] 085601材料工程315分求调剂 +3 yang_0104 2026-03-15 3/150 2026-03-15 10:58 by peike
[考研] 材料与化工(0856)304求B区调剂 +6 邱gl 2026-03-12 7/350 2026-03-13 23:24 by 邱gl
[考研] 281求调剂 +9 Koxui 2026-03-12 11/550 2026-03-13 20:50 by Koxui
[考研] 考研调剂 +4 芬达46 2026-03-12 4/200 2026-03-13 16:04 by ruiyingmiao
[考研] 材料301分求调剂 +5 Liyouyumairs 2026-03-12 5/250 2026-03-13 14:42 by JourneyLucky
信息提示
请填处理意见