Originally posted by hakekill at 2010-09-09 09:36:04:
你的子程序里面do 1030程序段里面使用了啊，算dis2的时候。

Originally posted by snoopyzhao at 2010-09-09 11:43:48:
你需要在 common 中指定 xc, yc, zc 数组的维数。

另外，你主程序中的 a(Tatom) 没有定义。

Originally posted by zyj8119 at 2010-09-09 11:53:25:

怎么定义？

Originally posted by snoopyzhao at 2010-09-09 15:21:01:


你问我吗？呵呵……

仔细看了一下程序，很多地方都有问题，呵呵……

首先 a 这个数组没有定义，从前后文看，a 似乎就是 name，但不太确定。

其次，xc, yc, zc 从 2.car 中赋值后，又被全部重新赋值为 ...

Originally posted by snoopyzhao at 2010-09-09 15:21:01:


你问我吗？呵呵……

仔细看了一下程序，很多地方都有问题，呵呵……

首先 a 这个数组没有定义，从前后文看，a 似乎就是 name，但不太确定。

其次，xc, yc, zc 从 2.car 中赋值后，又被全部重新赋值为 ...

Originally posted by zyj8119 at 2010-09-13 09:54:56:

Tatom的最大值不是num吧？

      PROGRAM MAIN
      PARAMETER(num=20000,blx=61.634002,bly=62.406401,blz=63.429759,
     *        namino=200)
        CHARACTER name(num)*9,fft(num)*5,atomname(num)*4,fx(num)*4
        INTEGER occupation(num),Tatom,itt,i
      dimension xc,yc,zc,charge
         real xc(num),yc(num),zc(num),charge(num)
      real xtemp,ytemp,ztemp,xfinal,yfinal,zfinal
      OPEN(10,file='2.car',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)
        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)
             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)
                  
      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)
                  
      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)
                  
      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)
       
        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 j=1,14
        lct=sqrt((xtemp-center(i,1))**2+(ytemp-center(i,2)**2))
        if(lct.lt.lc)then
        lc=lct
*******把整个体分解成无数的平行的平面，确定整个平面的每个分割区域的中心*******************      
        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
*****************向负方向搜索**********************************
      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      
      

Originally posted by zyj8119 at 2010-09-13 09:57:36:

我的这个程序就是，在一个三维晶胞中生成随机位置，只要这个位置与晶胞中的原子的距离小于一定的数值就可以认为成键了，
[code]      PROGRAM MAIN
      PARAMETER(num=20000,blx=61.634002,bly=62.406401,b ...