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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » Monte Carlo » 【讨论】我收集到的状态方程编程
81/1返回列表
查看: 2253  |  回复: 7
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

zyj8119

木虫 (著名写手)

[交流] 【讨论】我收集到的状态方程编程 已有1人参与

CODE:
// hardspherechain.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include "iostream.h"
#include "fstream.h"
#include "stdlib.h"


#define N 240 //将狭缝分为N个区间
#define DIA 20//直径占的区间数
#define STEP 1800 //步数
#define YITA  0.05//对比密度
#define PI 3.141592653589796234
#define T 1.
#define L 1.5   //ramit
#define M 1 //链长
#define BV 0//链长
#define fix 0//从0开始到M-1
#define DELTA 0.00000001  //误差
#define K 250//链长
#define LL 1
#define P1 (-890.366+2510.86*sqrt(L)-2629.19*L+1212.91*L*sqrt(L)-208.167*L*L)
#define P2 (-957.906+2722.24*sqrt(L)-2872.35*L+1334.96*L*sqrt(L)-230.768*L*L)
#define P3 (943.572-2808.87*sqrt(L)+3082.31*L-1481.70*L*sqrt(L)+263.780*L*L)

#define deP1 (2510.86/2./sqrt(L)-2629.19+1.5*1212.91*sqrt(L)-2.*208.167*L)
#define deP2 (2722.24/2./sqrt(L)-2872.35+1.5*1334.96*sqrt(L)-2.*230.768*L)
#define deP3 (-2808.87/2./sqrt(L)+3082.31-1.5*1481.70*sqrt(L)+2.*263.780*L)


#define PA ((P1+P2+6.)/4.)
#define PB ((P1-P2+P3+5.)/7.)
#define PC ((2.*P1-P3+2.)/3.)


#define W(z) (0.75/DIA/DIA/DIA*(DIA*DIA-(z)*(z)))   //硬球权重因子
#define W2(z) (0.75/DIA/DIA/DIA/(L)/(L)/(L)*(DIA*DIA*(L)*(L)-(z)*(z)))//方井权重因子

#define A (M*(1+(M-1)*0.45696/(M*1.)-(M-1)*(M-2)*0.74745/(M*M*1.)))
#define B (M*(1+(M-1)*2.10386/(M*1.)+(M-1)*(M-2)*3.49695/(M*M*1.)))
#define C (M*(1+(M-1)*1.75503/(M*1.)+(M-1)*(M-2)*4.83207/(M*M*1.)))

#define BA(yita)  ((((3.+A-B+3.*C)*(yita)-(1.+A+B-C))/(2.*(1.-(yita)))+(1.+A+B-C)/(2*(1.-(yita))*(1.-(yita)))+(C-1.)*log(1.-(yita)))/M)
#define DEBA(yita) ((PI/6.*((1.-C)/(1.-(yita))+(1.-B+2.*C)/(1.-(yita))/(1.-(yita))+(1.+A+B-C)/(1.-(yita))/(1.-(yita))/(1.-(yita))))/M)

#define METHOD  0  //计算方法，0表示单链模拟，1表示GL、GR方法

double rou[N],newrou[N],rouav[N],rouav2[N],su[M];   
double lamda[N],lamda2[N],gl[M][N],gr[M][N];
double famda[M][N],fvext[M][N],roum[M][N];
double avgdens,correct;   //平均密度
double dd;  //误差


double monoba(double yita)
{
        double khs,a1,a2,i1,i2;
        khs=(1.-yita)*(1.-yita)*(1.-yita)*(1.-yita)/(1.+4.*yita+4.*yita*yita-4.*yita*yita*yita+yita*yita*yita*yita);
        /////////////
        i1=(P1*yita+P2)/2./yita/(1.-yita);
        i1+=-P2/2./yita/(1.-yita)/(1.-yita);
        i1+=P3*log(1.-yita)+1.;
        /////////////////////
        i2=P1/2./(1.-yita);
        i2+=(P1*yita+P2)/2./(1.-yita)/(1.-yita);
        i2+=-P2/(1.-yita)/(1.-yita)/(1.-yita);
        i2+=P3*log(1.-yita);
        i2+=-yita*P3/(1.-yita)+1.;
        //i2=1/3*(yita*yita*yita-1)*i2;
        /////////////////////
        a1=-4.*(L*L*L-1.)/T*yita*i1;
        a2=-2.*(L*L*L-1.)/T/T*yita*khs*i2;
        return (a1+a2);
}


double demonoba(double yita)     //方阱单体项的导数
{
        double i2,khs;
        double da1dy,da2dy,dkhsdy,di2dy,dyitai2dy;
        khs=(1.-yita)*(1.-yita)*(1.-yita)*(1.-yita)/(1.+4.*yita+4.*yita*yita-4.*yita*yita*yita+yita*yita*yita*yita);
        /////////////
        /////////////////////
        i2=P1/2./(1.-yita);
        i2+=(P1*yita+P2)/2./(1.-yita)/(1.-yita);
        i2+=-P2/(1.-yita)/(1.-yita)/(1.-yita);
        i2+=P3*log(1.-yita);
        i2+=-yita*P3/(1.-yita)+1.;
        /////////////////////
        dkhsdy=-4.*(1.-yita)*(1.-yita)*(1.-yita)/(1+4.*yita+4.*yita*yita-4.*yita*yita*yita+yita*yita*yita*yita)-(1.-yita)*(1.-yita)*(1.-yita)*(1.-yita)*(4.+8.*yita-12.*yita*yita+4.*yita*yita*yita)/(1+4.*yita+4.*yita*yita-4.*yita*yita*yita+yita*yita*yita*yita)/(1.+4.*yita+4.*yita*yita-4.*yita*yita*yita+yita*yita*yita*yita);
        di2dy=(P1-P3)/(1.-yita)/(1.-yita)+((1.-yita)*(P1*yita+P2)-3.*P2)/(1.-yita)/(1.-yita)/(1.-yita)/(1.-yita)-P3/(1.-yita);
        dyitai2dy=i2+yita*di2dy;
        da1dy=-4.*(L*L*L-1.)/T*i2;
        da2dy=-2.*(L*L*L-1.)/T/T*(dkhsdy*yita*i2+khs*dyitai2dy);
        return ((da1dy+da2dy)*PI/6.);
}






double chainba(double yita)
{
        double ghs,i1dl,i1,i2,a1dl,a1dy,lny1,lny2;
        ghs=(1.-1./2.*yita)/(1.-yita)/(1.-yita)/(1.-yita);
        /////////////
        i1=(P1*yita+P2)/2./yita/(1.-yita);
        i1+=-P2/2./yita/(1.-yita)/(1.-yita);
        i1+=P3*log(1.-yita)+1.;
                /////////////////////       
        i1dl=(deP1*yita+deP2)/2./yita/(1.-yita);
        i1dl+=-deP2/2./yita/(1.-yita)/(1.-yita);
        i1dl+=deP3*log(1.-yita);
        /////////////////////
        i2=P1/2./(1.-yita);
        i2+=(P1*yita+P2)/2./(1.-yita)/(1.-yita);
        i2+=-P2/(1.-yita)/(1.-yita)/(1.-yita);
        i2+=P3*log(1.-yita);
        i2+=-yita*P3/(1.-yita)+1.;
        ////////////////////////////////
        a1dy=-4.*(L*L*L-1.)/T*i2;
        a1dl=-4.*(L*L*L-1.)/T*yita*i1dl-4.*3.*L*L/T*yita*i1;
        /////////////////////

lny1=log((ghs+1./4*(a1dy-L/3./yita*a1dl))/ghs);
//////////////////////////////////////

        lny2=(M-1.)/M/T*((PA*yita+PB)/2./(1.-yita)-PB/2./(1.-yita)/(1.-yita)+PC*log(1.-yita));

        return (-(M-1.)*lny1-(M-2.)*lny2);
}



double dechainba(double yita)
{
        return (chainba(yita+0.0001)-chainba(yita))/0.0001*PI/6.;
}

///////////////////////////////////////////////////////////////////////////
double ba(double yita)   
{
        return (monoba(yita)+chainba(yita)/M);
}


double deba(double yita)   
{
        return (demonoba(yita)+dechainba(yita)/M);
}

int min(int i,int j)
{
        if(i>j)return j;
        else return i;
}

int max(int i,int j)
{
        if(i>j)return i;
        else return j;
}

double mind(double i,double j)
{
        if(i         else return j;
}

double maxd(double i,double j)
{
        if(i>j)return i;
        else return j;
}

double juedui(double x)
{
        if(x<0)return -x;
        return x;
}

double bvext(int z)
{
        if(z=N-DIA)return 1.E+100;
        return 0;
        //return 0;
}
double bfix(int z)
{
        if(z=N-DIA)return 1.E+100;
        if(z<(DIA+LL))return -BV;//方井壁
        //return -15*(exp(-LAM*(z-DIA)*1./DIA)+exp(-LAM*(N-DIA-z)*1./DIA));//指数壁
        //if(z<(DIA+LL*DIA)) return BV*z/(LL*DIA)-(LL+1)*BV/LL;//三行一起用 表示三角壁
        //if(z>=(N-DIA-LL*DIA)) return BV*(N-z)/(LL*DIA)-(LL+1)*BV/LL;
        return 0;
}
void firstrou()     //初始化密度分布
{
        int i;
        double sum;
        for(i=0,sum=0.;i         {
                if(bvext(i)>1.E+90)rou[i]=0.;
                else rou[i]=YITA*6./PI;
                sum+=rou[i];
        }
        avgdens=sum/(N-2.*DIA);
        return;
}

void calrouav()     //计算粗粒化密度
{
        int i,j,j1,j2;
        double sum;

        for(i=0;i         {
                sum=0.;
             j1=max(DIA,i-DIA);
                j2=min(N-1-DIA,i+DIA);
                sum+=W(i-j1)*rou[j1]/2.;
                for(j=j1+1;j                 {
                        sum+=W(i-j)*rou[j];
                }
                sum+=W(i-j2)*rou[j2]/2.;
               
                rouav[i]=sum;//
        }
        //////////////以下为修正项///////////////////
        //for(i=DIA,sum=0.;i         //correct=avgdens*(N-2.*DIA)/sum/2.;
        //for(i=DIA,sum=0.;i     /////////////////////////////////////////////
        for(i=0;i         {
                sum=0.;
                        j1=max(DIA,i-int(L*DIA));
                j2=min(N-1-DIA,i+int(L*DIA));

                sum+=W2(i-j1)*rou[j1]/2.;
                for(j=j1+1;j                 {
                        sum+=W2(i-j)*rou[j];
                }
                sum+=W2(i-j2)*rou[j2]/2.;
               
                rouav2[i]=sum;//
        }
        ////////////以下为修正项////////////////////
        //for(i=DIA,sum=0.;i         //correct2=avgdens*(N-2.*DIA)/sum/2.;
        //for(i=DIA,sum=0.;i         /////////////////////////////////////////////
        return;
}

void callamda()     //计算lamda
{
        int i,j,k,j1,j2;
        double sum;

        for(i=0;i         {
                sum=0.;
              j1=max(DIA,i-DIA);
                j2=min(N-1-DIA,i+DIA);
        sum+=rou[j1]*W(j1-i)*DEBA(rouav[j1]*PI/6.)/2.;
                for(j=j1+1;j                 sum+=rou[j2]*W(j2-i)*DEBA(rouav[j2]*PI/6.)/2.;
                //sum*=correct;//
                sum+=BA(rouav[i]*PI/6.);
                lamda[i]=sum;
        }
        for(i=0;i         {
                sum=0.;
                        j1=max(DIA,i-int(L*DIA));
                j2=min(N-1-DIA,i+int(L*DIA));
        sum+=rou[j1]*W2(j1-i)*deba(rouav2[j1]*PI/6.)/2.;
                for(j=j1+1;j                 sum+=rou[j2]*W2(j2-i)*deba(rouav2[j2]*PI/6.)/2.;
                sum+=ba(rouav2[i]*PI/6.);
                lamda2[i]=sum;//cout<         }
        if(METHOD==1)
        {
                for(i=0;i                 {  lamda[i]+=lamda2[i]+bvext(i);
                     for(k=0;k                          {fvext[k][i]=0.;
                        if(k==fix)fvext[k][i]=bfix(i);
                         famda[k][i]=lamda[i]+fvext[k][i];
                         }
                }
                for(i=0;i                 {
                        gl[0][i]=1.;
                        gr[M-1][i]=1.;
                }
                for(i=0;i                 {
           j1=max(DIA,i-DIA);
                   j2=min(N-1-DIA,i+DIA);
                   for(k=1;k                    {
                           gl[k][i]=exp(-famda[k-1][j1])*gl[k-1][j1]/2.;
                           gr[M-k-1][i]=exp(-famda[M-k][j1])*gr[M-k][j1]/2.;
                           for(j=j1+1;j                            {
                                   gl[k][i]+=exp(-famda[k-1][j])*gl[k-1][j];
                               gr[M-k-1][i]+=exp(-famda[M-k][j])*gr[M-k][j];
                           }
                           gl[k][i]+=exp(-famda[k-1][j2])*gl[k-1][j2]/2.;
                           gr[M-k-1][i]+=exp(-famda[M-k][j2])*gr[M-k][j2]/2.;
                           gl[k][i]=gl[k][i]/(2*DIA);
                           gr[M-k-1][i]=gr[M-k-1][i]/(2*DIA);
                   }
                }
        }
        if(METHOD==0)
        {
                for(i=0;i                 {  lamda[i]+=lamda2[i]+bvext(i);
                     for(k=0;k                          {fvext[k][i]=0.;
                        if(k==fix)fvext[k][i]=bfix(i);
                         famda[k][i]=lamda[i]+fvext[k][i];
                         }
                }

        }
        return;
}

void calrou()       //计算密度
{
        int i,z;
        if(METHOD==1)
        {
                for(z=0;z                 {
                        for(i=0;i                         {
                                roum[i][z]=exp(-famda[i][z])*gl[i][z]*gr[i][z];
                        }
                }
        }
                if(METHOD==0)
        {
                for(z=0;z                 {
                        for(i=0;i                         {
                                roum[i][z]=exp(-famda[i][z]);
                        }
                }
        }
        return;
}

int calnewrou()     //得到新的密度分布
{
        int i,k;
        double sum;
    for(k=0;k         { su[k]=0;
                for(i=0;i                 {       
                        su[k]+=roum[k][i];
                }
                su[k]=su[k]/(N-2.*DIA);
        }
        for(i=0;i         {
                for(k=0;k                 {
                        roum[k][i]=roum[k][i]*avgdens/su[k]/M;
                }
        }
        sum=0.;
        for(i=0;i         {sum=0;
            for(k=0;k                 {
            sum+=roum[k][i];
                }
                newrou[i]=sum;
        }

        dd=0.;
        for(i=0;i         {
                newrou[i]=(K*rou[i]+newrou[i])/(K+1.);
                dd=maxd(dd,juedui(newrou[i]-rou[i]));
                rou[i]=newrou[i];
        }
        if(dd         return 0;
}

void out()         //输出
{
        int i;
        ofstream Op("output.dat",ios::out);               
        for(i=DIA;i         {
                Op<         }
    cout<<"complite"<         return;       
}

int main(int argc, char* argv[])
{
        int i,condition;
        firstrou();
        for(i=0;i         {
                calrouav();
                callamda();
                calrou();
                condition=calnewrou();
                if(condition==1)break;
                if(i%100==0)cout<         }
        out();
        cout<         cout<         return 0;
}
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
好好学习，天天向上。
1楼 2010-09-08 17:58:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
★ ★ ★
lei0736(金币+3):谢谢 2010-09-08 19:20:48
引用回帖:
Originally posted by zyj8119 at 2010-09-08 17:58:47:
[code]// hardspherechain.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include "iostream.h"
#include &quo ...

其实，MC吸附也可以与链流体密度泛函的结果相互比对，再与实验的比对，这三者分别为模拟，理论，实验。在没有模拟前，只有理论和实验。
赞一下(1人) 回复此楼
好好学习，天天向上。
2楼2010-09-08 18:00:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
★ ★ ★
lei0736(金币+3):谢谢 2010-09-08 19:20:59
MOF或者分子筛的结构可以划分为无数的小的格子，使用状态方程，通过迭代，求出每个格子的最终吸附数量，然后加和，就是总的吸附数量。而模拟却只是考虑LJ势能函数和电荷，不必使用状态方程迭代而已，使用的是MC方法。
赞一下(1人) 回复此楼
好好学习，天天向上。
3楼2010-09-08 18:02:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
★ ★
lei0736(金币+2):谢谢 2010-09-08 19:21:09
对MOF或者分子筛吸附的流体密度泛函感兴趣，可以参考 胡英的《近代化工热力学》，以及姜建文的理论文章（非模拟）
赞一下(1人) 回复此楼
好好学习，天天向上。
4楼2010-09-08 18:02:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
★
lei0736(金币+1):谢谢 2010-09-08 19:21:19
CO2也是一种链流体，不过挺简单的。
赞一下(1人) 回复此楼
好好学习，天天向上。
5楼2010-09-08 18:05:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
★ ★ ★ ★ ★ ★
zh1987hs(金币+3):谢谢专家的资源，希望能够再接再厉 2010-09-10 23:09:11
zh1987hs(金币+3):谢谢专家的资源，希望能够再接再厉 2010-09-10 23:10:14
引用回帖:
Originally posted by zyj8119 at 2010-09-08 18:05:34:
CO2也是一种链流体，不过挺简单的。

http://d.namipan.com/d/ad63c9a09 ... a72998111dac15b0000：这个是双链铆接，上面的程序是方阱纯链铆接。
赞一下(2人) 回复此楼
好好学习，天天向上。
6楼2010-09-10 16:51:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
★ ★
zh1987hs(金币+1):谢谢专家的解释 2010-09-10 23:09:39
zh1987hs(金币+1):谢谢专家的解释 2010-09-10 23:10:04
引用回帖:
Originally posted by zyj8119 at 2010-09-08 17:58:47:
[code]// hardspherechain.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include "iostream.h"
#include &quo ...

其实做这个的也很少，我看到前面有个人要硬球的代码，如果能帮助到需要的人，我觉得就很好了，这些代码都是不同状态方程的近似解，要想搞懂，必须对近代化工热力学很熟悉的。。
赞一下(2人) 回复此楼
好好学习，天天向上。
7楼2010-09-10 18:13:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
★
zh1987hs(金币+1):谢谢 2010-09-10 23:09:52
引用回帖:
Originally posted by zyj8119 at 2010-09-08 17:58:47:
[code]// hardspherechain.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include "iostream.h"
#include &quo ...

我觉得能够看过代码的人，都可以七嘴八舌的说几句，说错了没有关系，只有这样大家的水平才能上升。
赞一下(1人) 回复此楼
好好学习，天天向上。
8楼2010-09-10 18:16:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zyj8119 的主题更新
81/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 今天务委会开完了，明天出结果吗 +17 angus9576 2026-08-25 21/1050 2026-08-25 22:39 by 悠悠美女
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件，有关系吗？ +12 木水思豆 2026-08-25 14/700 2026-08-25 22:35 by yys300
[基金申请] 在坚冰还盖着北海的时候，我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 明天应该可查了！？ +6 chengyan1220 2026-08-23 6/300 2026-08-25 19:45 by zfd97
[基金申请] 放榜前的不淡定 20+4 snowwithsea 2026-08-19 19/950 2026-08-25 17:57 by chick875
[基金申请] 2026国自然函评费到账 +20 羊腰板 2026-08-21 23/1150 2026-08-25 16:34 by zsna
[基金申请] 什么时候开奖？ +12 CrisMessi 2026-08-18 13/650 2026-08-25 14:40 by ZJTJZ
[基金申请] 如果此刻你正在为国基感到焦虑，不妨来听听这首《基金之外》 +8 scalable 2026-08-24 8/400 2026-08-25 12:52 by jnhyjjm
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
[基金申请] 人气不行了 +11 fansofjerry 2026-08-21 11/550 2026-08-25 11:04 by 孤独的英雄6
[教师之家] 导师吐槽：我怎么摊上了这么个极品研究生！ +3 苏东坡二世 2026-08-23 3/150 2026-08-25 10:35 by shisan1313
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +5 process2012 2026-08-23 7/350 2026-08-25 09:42 by huixian257
[基金申请] 今天基金会出结果吗？20260819 +17 kkkl_v 2026-08-19 18/900 2026-08-25 09:41 by windflowerwy
[基金申请] 只有每年这种时候来逛逛小木虫 +25 yaoyewhu2008 2026-08-20 27/1350 2026-08-25 08:01 by Equinoxhua
[基金申请] 能否退出参与的面上项目解除限项 +21 koalala 2026-08-24 24/1200 2026-08-24 19:25 by 家与远方
[基金申请] 建议基金发布提前给出明确的时间点 +13 kulium 2026-08-21 16/800 2026-08-24 16:27 by superceng
[基金申请] 科研孤儿太难了 +18 我4大白菜 2026-08-20 19/950 2026-08-24 09:47 by 凯恩广盛催化分�
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[基金申请] 应该是下周三26日公布了吧？ +4 哈哈蛤？ 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[基金申请] 今天放榜没戏了吧 +9 yuleib84 2026-08-19 11/550 2026-08-21 10:06 by gltch
信息提示
关闭
请填处理意见
关闭