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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » Monte Carlo » 【原创】MCM-41的孔径分布的计算
141/2返回列表12下一页
查看: 2856  |  回复: 13
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 1 个 模拟EPI ，点击这里进行查看

zyj8119

木虫 (著名写手)

[交流] 【原创】MCM-41的孔径分布的计算 已有4人参与

接上一个帖子，MCM-41建好的模型会有孔径分布，这个在定义中是很复杂的，应该怎么计算呢，可以导出MCM-41中所有的H原子，然后自己写程序去计算孔径分布，我使用的是FORTRAN语言，导出的car结构中，可以使用EXCEL找到所有的H原子，这些H原子的 car文件可以保存为一个名叫H.car的文件，这个文件就可以进行孔径分布的计算，car文件的优点是带有电荷。读入这个car文件，可以输出一个名叫output.car的文件。
CODE:
        integer hno
        parameter(hno=2356)

        character a(hno)*5,fx(hno)*4,fft(hno)*5,natom(hno)*1
        double precision xc(hno),yc(hno),zc(hno)
        integer occupation(hno)
        real charge(hno)


      open(10,file='H.car',status='old')
        do 20 i=1,hno
        read(10,*)a(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i),
     &          natom(i),charge(i)
20    continue
      close(10)

        jmm=0
      open(30,file='output.car',access='append')
*40        amm=RAN2(IDUM)
40        imm=int(RAN2(IDUM)*hno)
        if(charge(imm).lt.1.0)then
        jmm=jmm+1
        write(30,888)a(imm),xc(imm),yc(imm),zc(imm),fx(imm),
     &             occupation(imm),fft(imm),natom(imm),charge(imm)
        charge(imm)=2.0
        endif



        if(jmm.eq.50.or.jmm.eq.100.or.jmm.eq.150.or.jmm.eq.200.or.
     &          jmm.eq.250.or.jmm.eq.300.or.jmm.eq.350.or.jmm.eq.400)then
        write(30,*)'********************',jmm,'**************************'
        write(30,*)'********************',jmm,'**************************'
        endif

        if(jmm.lt.415)then
      goto 40
        else
        goto 50
        endif
50    close(30)

888   format(A5,3X,F12.9,2X,F13.9,3X,F12.9,1X,A4,1X,I1,6X,A5,3X,A1,F8.3)   
        end
       
      
      
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
*     Random Number Generator               
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

                          FUNCTION RAN2(IDUM)
  
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL ran2,AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
        INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END

其中的nho是你的H.car中H原子的个数，可以按照你的文件的不同随意更改成你所需要的值，得到的output.car可以做下一步的计算。

[ Last edited by zyj8119 on 2010-12-2 at 14:14 ]
回复此楼

» 猜你喜欢

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

qphll

金虫 (正式写手)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
赞一下.

自己动手, 其乐无穷.
赞一下(1人) 回复此楼
Life, Love, Laugh.
2楼2010-12-02 00:27:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
引用回帖:
Originally posted by qphll at 2010-12-02 00:27:09:
赞一下.

自己动手, 其乐无穷.

多谢，这个输出的output.car文件在下一步还有用，程序明天给出。
赞一下 回复此楼
好好学习，天天向上。
3楼2010-12-02 00:28:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
引用回帖:
Originally posted by zyj8119 at 2010-12-02 00:28:52:

多谢，这个输出的output.car文件在下一步还有用，程序明天给出。

嗯, 期待!

我也可以喘口气, 更新一下前面的几个帖子. 连着看了三天的lammps manual.....
赞一下(1人) 回复此楼
Life, Love, Laugh.
4楼2010-12-02 01:12:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
引用回帖:
Originally posted by zyj8119 at 2010-12-02 00:24:27:
接上一个帖子，MCM-41建好的模型会有孔径分布，这个在定义中式很复杂的，应该怎么计算呢，可以导出MCM-41中所有的H原子，然后自己写程序去计算孔径分布，我使用的是FORTRAN语言，导出的car结构中，可以使用EXCEL找 ...

今天开始进行下一步的工作，昨天得到的output.car是修正后的H的结构文件，如果要计算孔径分布，首先得统计孔径分布。统计孔径分布的程序：
CODE:
          integer noh
         open(10,file='noh.txt',status='old')
         read(10,*)noh
         close(10)
         call distribution(noh)
         end

         subroutine distribution(noh)
         integer noh,nr
         real a(noh,2),step,min,max,rt,rp,b(-4:44,2)
         step=0.5
         min=0.0
         max=0.0
         rt=0.0
         rp=0.0
         nr=0.0
         open(20,file='input.txt',status='old')
         do 30 i=1,noh
         read(20,*)(a(i,j),j=1,2)
30     continue
        close(20)

        do 40 j=-4,44
        min=j*step-0.25
        max=j*step+0.25
               do 50 i=1,noh
        if(min.le.a(i,1).and.max.gt.a(i,1))then
        nr=nr+1
        rt=rt+a(i,2)
        end if
50    continue
       rp=rt/(1.0*nr)
       b(j,1)=nr
       b(j,2)=rp
       nr=0
       rt=0
       rp=0
40   continue
      open(60,file='output.txt',access='append')
     write(60,*)'*****************amount of H in the region**************'
       do 70 i=-4,44
       write(60,*)b(i,1)
70   continue
       write(60,*)'*****************average radius of the region***************'
       do 80 i=-4,44
       write(60,*)b(i,2)
80   continue
      end subroutine

此处的noh也是H的原子个数，如果你不想打开noh.txt,直接给noh一个parameter语句也是可以的。

[ Last edited by zyj8119 on 2010-12-2 at 09:41 ]
赞一下 回复此楼
好好学习，天天向上。
5楼2010-12-02 09:20:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hsqlxsz

新虫 (著名写手)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
眼前一片茫然，技术含量太高，仰视一下
赞一下(1人) 回复此楼
6楼2010-12-02 09:34:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
统计完成以后，接下来就是计算孔径分布了，源程序是这样的：
CODE:
*               main program
                real cx,cy
               integer HNO,lattice
               open(1000,file='parameter.txt',status='old')
               read(1000,*)cx
                     read(1000,*)cy
                     read(1000,*)HNO
               read(1000,*)lattic
               call statistic(cx,cy,HNO,lattic)
*              write(*,*)cx,cy,HNO
               end

               subroutine statistic(cx,cy,HNO,lattic)
               real cx,cy,cz,fc0,cx0,cy0,range,step
               integer HNO,cyc,lattic
               real p(HNO,11),r(lattic,lattic,HNO),rp(lattic,100)
               real d(HNO),rmin,rmax,fc(lattic,lattic)
               integer nb(0:99)
               range=10.0
              open(10,file='input.txt',status='old')
              do 10 i=1,HNO
              read(10,*)(p(i,j),j=1,11)
*            if(p(i,6).gt.29)then
*            p(i,6)=p(i,6)-57.975988
*            end if
*            if(p(i,7).gt.67)then
*            p(i,7)=p(i,7)-133.758958
*            end if
100       continue
            close(10)
***************print the coordinate***********************
            open(1000,file='output.txt',access='append')
            write(1000,*)'**********x coordinate*****************
            do 1001 i=1,HNO
              write(1000,*)p(i,6)
1001     continue
            write(1000,*)'**********y coordinate*****************
               do 1002 i=1,HNO
              write(1000,*)p(i,7)
1002     continue
             write(1000,*)'**********y coordinate*****************
               do 1003 i=1,HNO
              write(1000,*)p(i,8)
1003     continue
            close(1000)
********************************************************
             do 220 cyc=1,10
             fc0=1.0e20
             cx0=cx
             cy0=cy
             step=range*2.0/(1.0*lattic)
             do 101 i=1,lattic
                do 102 j=1,lattic
              rp(i,j)=0.0
              fc(i,j)=0.0
102        continue
101        continue

             do 110 i=1,lattic
               do 120 j=1,lattic
                 do 130 k=1,HNO
                  cx=cx0-range+step*i
                  cy=cy0-range+step*j
                  r(i,j,k)=sqrt((p(k,6)-cx)**2+(p(k,6)-cy)**2)
                  rp(i,j)=rp(i,j)+r(i,j,k)/(HNO*1.0)
130            continue
120           continue
110         continue

              do 140 i=1,lattic
               do 150 j=1,lattic
                 do 160 k=1,HNO
                fc(i,j)=fc(i,j)+(r(i,j,k)-rp(i,j))**2
160           continue
               fc(i,j)=sqrt(fc(i,j)/1.0*(HNO-1)))
               if(fc(i,j).lt.fc0)then
               fc0=fc(i,j)
               cx=cx0-range+step*i
               cy=cy0-range+step*j
               end if
150         continue
140         continue
              range=range/2.0
220        continue
            
              do 170 i=1,HNO
                d(i)=sqrt((p(i,6)-cx)**2+(p(i,7)-cy)**2)
170         continue

              open(20,file='output.txt',access='append')
              write(20,*)'the center of pore'
              write(20,*)'x=',cx,'y=',cy
              write(20,*)'fangcha=',fc0
              write(20,*)'****************************'
              write(20,*)'distance between Hydrogen atom and the center of pore'
              do 180 i=1,HNO
              write(20,*)d(i)
180        continue
              close(20)
              do 185 i=0,99
              nb(i)=0
185        continue
     
             do 190 i=0,99
             rmin=10.0+0.1*i
             rmax=10.0+0.1*(i+1)
                 do 200 j=1,HNO
                   if(rmin.lt.d(j).and.rmax.ge.d(j))then
                  nb(i)=nb(i)+1
200          continue
190          continue

             open(30,file='output.txt',access='append')
             write(30,*)'****************nb*****************'
             write(30,*)'distribution of pore radius'
             do 210 i=0,99
             rr=10.1+0.1*i
             write(30,*)nb(i)
*           write(30,*)nb(i)
210        continue
             close(30)
             end
           

   

对文件中nb(i)画图拟合，就是所要得到的孔径分布曲线。

[ Last edited by zyj8119 on 2010-12-2 at 10:03 ]
赞一下 回复此楼
好好学习，天天向上。
7楼2010-12-02 10:01:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
引用回帖:
Originally posted by zyj8119 at 2010-12-02 10:01:19:
统计完成以后，接下来就是计算孔径分布了，源程序是这样的：

[code]
*               main program
                real cx,cy
               integer HNO,lattice
               open(1000,file='para ...

孔径分布都是从三维坐标中找到二维的进行计算。
赞一下 回复此楼
好好学习，天天向上。
8楼2010-12-02 10:50:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
引用回帖:
Originally posted by zyj8119 at 2010-12-02 10:01:19:
统计完成以后，接下来就是计算孔径分布了，源程序是这样的：

[code]
*               main program
                real cx,cy
               integer HNO,lattice
               open(1000,file='para ...

版主请评阅下。。。
回复此楼
好好学习，天天向上。
9楼2010-12-02 11:17:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
引用回帖:
Originally posted by hsqlxsz at 2010-12-02 09:34:40:
眼前一片茫然，技术含量太高，仰视一下

认真看, 使劲问. 越是能垒高的地方, 越是能很快进步.
赞一下(1人) 回复此楼
Life, Love, Laugh.
10楼2010-12-02 12:45:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zyj8119 的主题更新
141/2返回列表12下一页
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 投票:  有多少人是今天查系统知道结果的？ +15 爱看书的可乐 2026-08-26 17/850 2026-08-27 07:12 by xxniao123
[教师之家] 售SCI一区文章，我:8.O.551.O.5.4,科目全,可伽急 +3 LIbGuocjEEYw 2026-08-26 4/200 2026-08-27 04:33 by Ie9AyIAvGbvs
[基金申请] 看板上这么多中的，有点像50人群里49个人都是骗子的那种感觉…… +3 a089 2026-08-26 3/150 2026-08-26 23:21 by young90
[基金申请] 我不理解！ +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 怎么查啊 +5 huang1991js 2026-08-26 5/250 2026-08-26 16:09 by Equinoxhua
[基金申请] 范进中举一文的中心思想 +9 炎黄贵胄 2026-08-22 10/500 2026-08-26 15:40 by semaglutide
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件，有关系吗？ +25 木水思豆 2026-08-25 28/1400 2026-08-26 14:53 by draco1987
[基金申请] 哪位高人中了，把查询到的截图贴出来让我看看，让我长长见识 +4 yuleib84 2026-08-26 5/250 2026-08-26 13:40 by yuleib84
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 国际合作可查了，中了面上 +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 怎么看青基中了没有啊 +4 叶九微 2026-08-26 4/200 2026-08-26 10:52 by xiacongxu
[职场人生] 学生撞见辅导员送外卖，第二天全班都沉默了 +3 鲸鱼融金_浙江_� 2026-08-22 3/150 2026-08-26 09:01 by zzuzxg
[基金申请] 牛来！米来！面来！ +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
[基金申请] 明天应该可查了！？ +6 chengyan1220 2026-08-23 6/300 2026-08-25 19:45 by zfd97
[基金申请] 某些机构，以效率低为荣，以效率低作为存在感 +9 yuleib84 2026-08-25 10/500 2026-08-25 17:14 by alexon
[基金申请] 今日不放榜？网传国自然预计 8 月 27 日可查结果 +17 医学老男孩 2026-08-20 22/1100 2026-08-25 15:36 by 医学老男孩
[基金申请] 让我中一个面上吧！ +13 大萍1987 2026-08-20 16/800 2026-08-24 10:23 by 太傻了
[基金申请] 看来今天不会放榜了？ +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
信息提示
关闭
请填处理意见
关闭