应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 其他 » 最优化计算求助！（1stopt或者Matlab）
51/1返回列表
查看: 1575  |  回复: 16
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

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 ]
回复此楼

» 猜你喜欢

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

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2012-03-26 13:54:34
已阅   回复此楼   关注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的回帖
查看全部 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角度后（第三个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 个回答
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭