program main
    parameter(natom=12348,nN=1000)
    integer atomsn,num(natom)
    real xc,yc,zc,dl,coor(natom,3),zuo(nN,3)
    character atomlabel(natom)*2,tlabel(natom)*2
      open(10,file='1.txt',status='old')
    do 15 j=1,natom
    read(10,*)num(j),atomlabel(j),coor(j,1),coor(j,2),coo
     &    r(j,3),tlabel(j)
15    continue
      close(10)
C     generate the initiated coordinates for nitrogen atom
C     check the distance between nitrogen and atoms in framework
      do 25 i=1,nN
     do 26 j=1,3
    zuo(i,j)=0
26    continue
25    continue

      do 30 i=1,nN
50    xc=RAN2(IDUM)*58.015
      yc=RAN2(IDUM)*133.80
    zc=RAN2(IDUM)*39.168
       do 40 j=1,natom
       dl=sqrt((xc-coor(j,1))**2+(yc-coor(j,2))**2+(zc-coor(j,3))**2)
       if(dl.lt.1.9)goto 50
40    continue
C     check the distance between nitrogen atoms
      if(i.ge.2)then
    do 60 k=1,i-1
     dl=sqrt((xc-zuo(k,1))**2+(yc-zuo(k,2))**2+(zc-zuo(k,3))**2)
     if(dl.lt.1.9)goto 50
60    continue
      end if
      zuo(i,1)=xc
    zuo(i,2)=yc
    zuo(i,3)=zc
30    continue
C     ouput the file for next simulation      
    open(70,file='2.txt',access='append')
      NAME='N'
C     output the coordinates of nitrogen    
    do 120 j=1,natom
        write(70,666)num(j),atomlabel(j),coor(j,1),coor(j,2),
     & coor(j,3),tlabel(j)
120   continue
C     output the coordinates of atoms in framework      
    do 140 i=1,nN
       atomsn=natom+i
      write(70,666)atomsn,NAME,zuo(atomsn,1),zuo(a
     & tomsn,2),zuo(atomsn,3),NAME
140   continue

666   format(I5,A2,3X,F9.3,2F7.3,3X,A2)
      close(70)

    end


C     obtain the randam number between 0 and 1    
    FUNCTION RAN2(INUM)
    INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
    REAL ran2,AM,eps,RNMX
    PARAMETER(IM1=214783563,IM2=214783399,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)
    end if
    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
    数据文件格式是:
    1    SI            15.65    37.172    15.303            Si
2    SI            15.762    56.005    14.655            Si
3    SI            14.211    74.922    14.369            Si
4    SI            16.512    94.068    14.554            Si
5    SI            14.894    131.378    15.386            Si
6    SI            34.373    0.61    16.75            Si
7    SI            53.776    17.43    16.264            Si
8    SI            53.474    114.575    15.138            Si
9    SI            11.633    4.363    14.668            Si
10    SI            14.412    24.617    13.956            Si
11    SI            13.28    43.362    13.062            Si
12    SI            13.565    82.058    12.901            Si
13    SI            12.393    121.022    12.86            Si
14    SI            33.626    102.581    12.622            Si
        运行不出来。。。