[½»Á÷] ¡¾ÇóÖú¡¿FORTRAN±àÒëͨ¹ý£¬µ«ÊÇÔËÐв»ÁË ÒÑÓÐ3È˲ÎÓë

PROGRAM MAIN
      PARAMETER(num=14797,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
      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)
      do 30 i=1,num
          xc(i)=0
          yc(i)=0
          zc(i)=0
30    continue
     
      natom_add=0
        Tatom=num+natom_add
*********************************************************
*     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)a(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)a(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)a(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)a(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)a(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)a(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
        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-08 16:51:55:
center.txt:
2.201  2.228
3.301  3.342
5.502  5.570
7.703  7.798
9.904  10.026
12.105 12.254
14.306 14.482
16.507 16.710
18.708 18.938
20.909 21.166
23.110 23.394
25.311 25.622
27.513 2 ...