[交流] 【求助】老虎大王的RDF计算程序已有2人参与

Originally posted by lei0736 at 2010-09-27 13:44:17:
您的文件已经成功上传
文件下载地址:http://g.zhubajie.com/urllink.php?id=9971798iop5bn34l37l5kf4
您可以通过复制和粘贴以下的链接地址，与您的朋友共享您的文件。

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!构型文件格式要求(设总原子数为ntotal，总共有nconfig个构型)
!!  A.第一行： a   b   c   (MD元胞边长)
!!  B. 第二行至第ntotal+1行： ind  x  y  z  (ind为原子种类，用1,2...表示), x范围从0到a, (y,z相同)
!!  C. 将B重复nconfig次
!! 如不符此格式，请自行修改读入语句
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!此行设置要根据情况修改!!!!!!!!!!!!!!
! nspec 是原子种类数, n_config是存贮的构型数

        parameter(nspec=1,n_config=1000)              

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        character*50 fname, filerdf
        dimension n(nspec)
        integer,dimension(:,:),allocatable:: nrdf
!        real,dimension(:,:),allocatable:: grdf       ! 注意grdf这里用不到，但如果要一起输出结果, 则很有用
        real,dimension(:),allocatable:: x
        real,dimension(:),allocatable:: y
        real,dimension(:),allocatable:: z
        integer,dimension(:),allocatable::ind

        data pi/3.14159265/

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!以下设置要根据情况更改

        data n/864/                !数组n是各类原子数目，对多元系，应是如 data n/864,864/ 之类
      fname='./myconf.dat'         !fname 构型文件名
        filerdf='./myrdf.rdf'          !filerdf 输出RDF的文件名
       rdflim=10.                !rdflim  RDF计算的长度(不要超过边长的一半)
        divrdf=0.01                !Divrdf  RDF曲线间隔
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        ntotal=0
        do i=1,nspec
        ntotal=ntotal+n(i)
        enddo
        inter=nspec*(nspec+1)/2
      conrdf=1/divrdf
        rcutd=rdflim
        maxint=1+int(conrdf*rcutd+0.5)
        rcutsd=rcutd*rcutd       
               
        write(*,*) ntotal
        write(*,*) inter
       
        allocate(x(ntotal))
        allocate(y(ntotal))
        allocate(z(ntotal))
        allocate(ind(ntotal))
        allocate(nrdf(maxint,inter))
!        allocate(grdf(maxint,inter))

               
       
      open(1,file=fname,form='formatted',status='old')
        read(1,*) a,b,c
      a=a/2
        b=b/2
        c=c/2
        write(*,*) a,b,c

       
        k=0
        do i=1,nspec
        do j=i,nspec
        k=k+1
        do jdiv=1,maxint
        nrdf(jdiv,k)=0
        enddo
        enddo
        enddo

        nsrd=0
        do kq=1,n_config
        if(mod(kq,200).eq.0) write(*,*) kq
        nj=0

        do j=1,ntotal
        read(1,*) ind(j),x(j),y(j),z(j)
        x(j)=x(j)-a
        y(j)=y(j)-b
        z(j)=z(j)-c
        enddo

        nsrd=nsrd+1
        nopm1=ntotal-1
        do i=1,nopm1
        ispec=ind(i)
        icnn=nspec*(ispec-1)-(ispec**2-ispec)/2
        ip1=i+1
        do j=ip1,ntotal
        rxij=x(i)-x(j)
        ryij=y(i)-y(j)
        rzij=z(i)-z(j)
        rx=rxij-2*a*aint(rxij/a)
        ry=ryij-2*b*aint(ryij/b)
        rz=rzij-2*c*aint(rzij/c)
        rsq=rx*rx+ry*ry+rz*rz
        if(rsq.gt.rcutsd) cycle  
        r=sqrt(rsq)
        jspec=ind(j)
        icn=icnn+jspec
        ibox=1+int(conrdf*r+0.5)
        nrdf(ibox,icn)=nrdf(ibox,icn)+1
        enddo
        enddo
       
       
        enddo            !end kq


!   calculate RDF       
        open(10,file=filerdf,status='new',form='formatted')

        icn=0
        do is=1,nspec
        do js=is,nspec
        icn=icn+1
        write(10,'(/1x,a,1x,i2,1x,i2)') 'RDF (distance in A) ', is,js
        fncrdf=float(nsrd*n(is))
        facrdf=8*a*b*c/(float(n(js))*4.*pi)
        qint=0.0
        maxmin=maxint-1
       
        do jdiv=2,maxmin
        qr=divrdf*float(jdiv-1)
        fnrdf=float(nrdf(jdiv,icn))
        if(is.eq.js)fnrdf=2.0*fnrdf
        qrdf=facrdf*fnrdf/(fncrdf*qr*qr*divrdf)
        qint=qint+fnrdf/fncrdf
        write(10,'(f8.3,2(1x,f11.3))')qr,qrdf,qint       !输出结果 r, RDF, 配位数
  !        grdf(jdiv,icn)=qrdf
        enddo

        enddo
        enddo
        close(1)
        close(10)

      end