Originally posted by zyj8119 at 2010-09-13 14:45:13:

有做過的蟲子，一起討論下。

Originally posted by zyj8119 at 2010-09-13 14:25:00:
使用MS，已經(jīng)建好了MCM-41的模型，現(xiàn)在想在MCM-41的內(nèi)部加上一定數(shù)目的SiCH2CH2CH2NH2這樣的結(jié)構(gòu)，用的是有限元分析法，就是把整個晶胞看成是一個立方體，導(dǎo)出為car文件，然后切割成14*14*14的小立方體，利用搜索 ...

Originally posted by zyj8119 at 2010-09-13 14:44:31:

這個問題困擾我很久了，對于MFI/MCM-41的內(nèi)部胺基鉚接問題，也是同類似的做，但是程序一直編不好。

Originally posted by zyj8119 at 2010-09-13 14:48:03:

http://www.gaoyang168.com/bbs/viewthread.php?tid=2299954

      PROGRAM MAIN
      PARAMETER(num=14797,blx=61.634002,bly=62.406401,blz=63.429759,
     *        namino=200)
        CHARACTER name(num),fft(num),atomname(num),fx(num)
        INTEGER occupation(num),Tatom,itt,i
      real xc(num),yc(num),zc(num),charge(num)
      real xtemp,ytemp,ztemp,xfinal,yfinal,zfinal
      OPEN(10,file='2.txt',status='old')
        do 20 i=1,num
          read(10,*)name(i),xc(i),yc(i),zc(i),fx(i),occupation(i),
     *  fft(i),atomname(i),charge(num)
20    continue
      close(10)
     
     
     
*********************************************************
*     add Si to the chosen Oxygen atom
*********************************************************
      do 45 itt=1,namino
        xtemp=blx*RAN2(IDUM)
        ytemp=bly*RAN2(IDUM)
        ztemp=blz*RAN2(IDUM)
        lbond=1.90
      call addatom(xtemp,ytemp,ztemp,lbond,Tatom,xfinal,yfinal,zfinal,xc
     *        ,yc,zc,num)
        natom_add=natom_add+1
      Tatom=Tatom+natom_add
      name(Tatom)='Si'
        xc(Tatom)=xfinal
        yc(Tatom)=yfinal
        zc(Tatom)=zfinal
        fx(Tatom)='XXXX'
        occupation(Tatom)=1
        fft(Tatom)='Si3'
      atomname(Tatom)='Si'
        charge(Tatom)=5.000
        open(140,file='amino.car',access='append')
        write(140,888)name(Tatom),xc(Tatom),
     &              yc(Tatom),zc(Tatom),
     &              fx(Tatom),occupation(Tatom),
     &              fft(Tatom),atomname(Tatom),
     &              charge(Tatom)
        close(140)

******************************************************************
*          add C1 atom on the chosen Si atom
******************************************************************
      lbond=1.93
        xtemp=xc(Tatom)
        ytemp=yc(Tatom)
        ztemp=zc(Tatom)
      call addatom(xtemp,ytemp,ztemp,lbond,Tatom,xfinal,yfinal,zfinal,xc
     *        ,yc,zc,num)
             natom_add=natom_add+1
        Tatom=Tatom+natom_add
      name(Tatom)='C'
        xc(Tatom)=xfinal
      yc(Tatom)=yfinal
        zc(Tatom)=zfinal
        fx(Tatom)='XXXX'
        occupation(Tatom)=1
        fft(Tatom)='C_3'
      atomname(Tatom)='C'
        charge(Tatom)=6.000

        open(150,file='amino.car',access='append')
        write(150,888)name(Tatom),xc(Tatom),
     &              yc(Tatom),zc(Tatom),
     &              fx(Tatom),occupation(Tatom),
     &              fft(Tatom),atomname(Tatom),
     &              charge(Tatom)
        close(150)
******************************************************************
*          add C2 atom on the chosen C1 atom
******************************************************************
      lbond=1.50
        xtemp=xc(Tatom)
        ytemp=yc(Tatom)
        ztemp=zc(Tatom)
      call addatom(xtemp,ytemp,ztemp,lbond,Tatom,xfinal,yfinal,zfinal,xc
     *        ,yc,zc,num)
                  
      natom_add=natom_add+1
        Tatom=Tatom+natom_add
      name(Tatom)='C'
        xc(Tatom)=xfinal
        yc(Tatom)=yfinal
        zc(Tatom)=zfinal
        fx(Tatom)='XXXX'
        occupation(Tatom)=1
        fft(Tatom)='C_3'
      atomname(Tatom)='C'
        charge(Tatom)=7.000

        open(160,file='amino.car',access='append')
        write(160,888)name(Tatom),xc(Tatom),
     &              yc(Tatom),zc(Tatom),
     &              fx(Tatom),occupation(Tatom),
     &              fft(Tatom),atomname(Tatom),
     &              charge(Tatom)
        close(160)
******************************************************************
*          add C3 atom on the chosen C2 atom
******************************************************************
      lbond=1.50
        xtemp=xc(Tatom)
        ytemp=yc(Tatom)
        ztemp=zc(Tatom)
          call addatom(xtemp,ytemp,ztemp,lbond,Tatom,xfinal,yfinal,zfinal,xc
     *        ,yc,zc,num)
                  
      natom_add=natom_add+1
        Tatom=Tatom+natom_add
      name(Tatom)='C'
        xc(Tatom)=xfinal
        yc(Tatom)=yfinal
        zc(Tatom)=zfinal
        fx(Tatom)='XXXX'
        occupation(Tatom)=1
        fft(Tatom)='C_3'
      atomname(Tatom)='C'
        charge(Tatom)=8.000

        open(170,file='amino.car',access='append')
        write(170,888)name(Tatom),xc(Tatom),
     &              yc(Tatom),zc(Tatom),
     &              fx(Tatom),occupation(Tatom),
     &              fft(Tatom),atomname(Tatom),
     &              charge(Tatom)
        close(170)
******************************************************************
*          add N atom on the chosen C3 atom
******************************************************************

      lbond=1.53
        xtemp=xc(Tatom)
        ytemp=yc(Tatom)
        ztemp=zc(Tatom)
             
          call addatom(xtemp,ytemp,ztemp,lbond,Tatom,xfinal,yfinal,zfinal,
     &          xc,yc,zc,num)
                  
      natom_add=natom_add+1
        Tatom=Tatom+natom_add
      name(Tatom)='N'
        xc(Tatom)=xfinal
        yc(Tatom)=yfinal
        zc(Tatom)=zfinal
        fx(Tatom)='XXXX'
        occupation(Tatom)=1
        fft(Tatom)='N_3'
      atomname(Tatom)='N'
        charge(Tatom)=9.000

        open(180,file='amino.car',access='append')
        write(180,888)name(Tatom),xc(Tatom),
     &              yc(Tatom),zc(Tatom),
     &              fx(Tatom),occupation(Tatom),
     &              fft(Tatom),atomname(Tatom),
     &              charge(Tatom)
        close(180)
45    continue
*********************************************************************************
*     write the output file
*********************************************************************************

        open(190,file='output-final.car',access='append')
        do 200 i=1,Tatom
        write(190,888)name(i),xc(i),yc(i),zc(i),fx(i),occupation(i),fft(i),
     &            atomname(i),charge(i)
200   continue
      close(190)


888   format(A5,3X,F12.9,2X,F13.9,3X,F12.9,1X,A4,1X,I1,6X,A5,3X,A2,F7.3)

*********************************************************************************
*********************************************************************************
                                        
      end



**********************************add atom************************************
        SUBROUTINE addatom(xtemp,ytemp,ztemp,lbond,Tatom,xfinal,yfinal,
     *        zfinal,xc,yc,zc,num)
       
        REAL center(14,2)
        INTEGER Tatom,im
        open(35,file='center.txt',status='old')
        do 36 i=1,14
        read(35,*)center(i,1),center(i,2)
36    continue
      close(35)
        lc=1000
        xcenter=0
        ycenter=0
        do 40 i=1,14
        lct=sqrt((xtemp-center(i,1))**2+(ytemp-center(i,2))**2)
        if(lct.lt.lc)then
        lc=lct
*******把整個體分解成無數(shù)的平行的平面，確定整個平面的每個分割區(qū)域的中心*******************      
        xcenter=center(i,1)
        ycenter=center(i,2)   
        end if
40    continue
      if(ytemp.le.ycenter)then
        ja=20
        jb=0
        jc=-1
        jd=j
        else
        ja=-20
        jb=0
        jc=1
        jd=j
        endif
*****************向負(fù)方向搜索**********************************
      if(xtemp.le.xcenter-5)then
        ia=20
        ib=0
        ic=-1
      minlbond=lbond-0.1
        maxlbond=lbond+0.1
      search_step=0.1
        do 1000 i=ia,ib,ic
          do 1010 j=ja,jb,jc
            do 1020 k=-20,20
            targetx=xtemp+search_step*i
            targety=ytemp+search_step*j
            targetz=ztemp+search_step*k
            dis1=sqrt((targetx-xtemp)**2+
     &              (targety-ytemp)**2+
     &              (targetz-ztemp)**2)
       if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then
            do 1030 im=1,Tatom
                   dis2=sqrt((targetx-xc(im))**2+
     &                (targety-yc(im))**2+
     &                (targetz-zc(im))**2)
                    if(dis2.lt.3.0)then
                                goto 1020
                   end if
             
1030  continue
      goto 1040
        end if
1020  continue
1010  continue
1000  continue
1040  xfinal=targetx
      yfinal=targety
        zfinal=targetz
        end if
*****************向正方向搜索**********************************
      if(xtemp.ge.xcenter+5)then
        ia=20
        ib=-20
        ic=-1
      minlbond=lbond-0.1
        maxlbond=lbond+0.1
      search_step=0.1
        do 1050 i=ia,ib,ic
          do 1060 j=ja,jb,jc
            do 1070 k=-20,20
            targetx=xtemp+search_step*i
            targety=ytemp+search_step*j
            targetz=ztemp+search_step*k
            dis1=sqrt((targetx-xtemp)**2+
     &              (targety-ytemp)**2+
     &              (targetz-ztemp)**2)
       if(dis1.gt.minlbond.and.dis1.lt.maxlbond)then
            do 1080 im=1,Tatom
                   dis2=sqrt((targetx-xc(im))**2+
     &                (targety-yc(im))**2+
     &                (targetz-zc(im))**2)
                    if(dis2.lt.3.0)then
                                goto 1070                   
                    endif
1080  continue
      end if
        goto 1090
1070  continue
1060  continue
1050  continue
1090  xfinal=targetx
      yfinal=targety
        zfinal=targetz
        end if                                    
      end SUBROUTINE
******************************************************************************
        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