24小时热门版块排行榜

返回列表

xk6891

至尊木虫 (著名写手)

模拟EPI: 1
应助: 0 (幼儿园)
金币: 13735.8
散金: 379
红花: 9
帖子: 1514
在线: 546小时
虫号: 890342
注册: 2009-11-01
专业: 理论和计算化学

[交流] 关于Radial Distribution Funciton的总结已有3人参与

原帖：http://hi.baidu.com/xk6891/blog/item/8b248511151749daa5ef3f6a.html
看了一些关于计算RDF的东西，曾经自己的思考也是如此，后来不知怎么忘了。现在总结一下。
查看的内容有：
1.http://www.phy.cmich.edu/people/petkov/isaacs/phys/rdfs.html
2.http://www.iams.sinica.edu.tw/lab/jlli/thesis_andy/node14.html
3.http://en.wikipedia.org/wiki/Radial_distribution_function
4.http://imagejdocu.tudor.lu/doku.php?id=macro:radial_distribution_function
5.http://muchong.com/bbs/viewthread.php?tid=2395138
6.http://www.materialssimulation.com/node/185
7.计算材料学基础（张跃.北航出版社.2007.）5.2.4节，P132
The Radial Distribution Function, R.D.F. , g(r), also called pair distribution function or pair correlation function, is an important structural characteristic, therefore computed by I.S.A.A.C.S..

Figure 1: Space discretization for the evaluation of the radial distribution function.
Considering a homogeneous distribution of the atoms/molecules in space, the g(r) represents the probability to find an atom in a shell dr at the distance r of another atom chosen as a reference point [Fig. 1].
By dividing the physical space/model volume into shells dr [Fig. 1] it is possible to compute the number of atoms dn(r) at a distance between r and r + dr from a given atom:
$dn(r) = (N/V)*g(r) 4πr^2dr (1)$
where N represents the total number of atoms, V the model volume and where g(r) is the radial distribution function.

上面是一个RDF函数的图像，横坐标是以一个原子/分子为中心与另一个原子/分子的距离，但纵坐标的理解是一个几率值，
而非在对应位置处找到另一个原子/分子的个数。张跃书中讲的比较清楚，“径向分布函数 g(r)是距离一个原子为 r 时找到另一个原子的概率 ,g (r)是一个量纲为 1的量。”可以看出g(r)是一个几率值而非数值，不然理想气体的径向分布函数g(r)=1，这里的1如何理解？难道理解为，在任何距离处找到另一个理想气体原子/分子的数目一直为1，这显然是不对的，因为不同密度的气体，每个粒子周围找到另一个粒子的数目会根据气体密度而变化(见ref.4)，所以这个1是一个几率值，只表示找到几率另一个粒子几率的大小。
再抄一个例子（见ref.6），以资佐证(注意下划线)：
对于很小的距离(小于原子间距) ,g (r)为0,这是由于强烈的排斥作用造成的。第一个(也是最高的 )峰出现在 r ≈ 3.7(埃) ,g (r) ≈ 3。这意味着两个分子在这个距离的几率是理想气体距离几率的三倍。然后径向分布函数下降,在 r ≈ 5.4(埃) 时达到最小值点 ,可知两原子在这个距离的几率比在理想气体状态时要小。随着距离的增加,g (r)趋近于理想气体时的值 1 ,这意味着体系不具有长程有序。
具体计算RDF的code已经写了很多，下面将收集到的几个例举如下(没有验证)：
java(参见ref.4):
*************************************************************************************

CODE:

BEGIN{

#Set a pair of atoms you want to consider.

     atoma="O";

     atomb="O";

#Set the size of cell: x, y, z.

     cellx=11.3061;

     celly=11.3061;

     cellz=11.3061;

#Set the minimum and maxinum dynamics steps you want to calculated.

     stepmin=50000;

     stepmax=60000;

#Set the interval (delta r, in Angstrom) and the total steps of radius (r).

     deltar=0.05;

     ntotal=200;

# Read user's initial parameter from the file "pair.txt" in current directory.

     for (i=1;i<=100;i++)

       {

       getline mypar < "pair.txt";

       split(mypar,mypar0);

       if ( mypar0[1] == "atoma" ) atoma=mypar0[2];

       if ( mypar0[1] == "atomb" ) atomb=mypar0[2];

       if ( mypar0[1] == "cellx" ) cellx=mypar0[2];

       if ( mypar0[1] == "celly" ) celly=mypar0[2];

       if ( mypar0[1] == "cellz" ) cellz=mypar0[2];

       if ( mypar0[1] == "stepmin" ) stepmin=mypar0[2];

       if ( mypar0[1] == "stepmax" ) stepmax=mypar0[2];

       if ( mypar0[1] == "deltar" ) deltar=mypar0[2];

       if ( mypar0[1] == "ntotal" ) ntotal=mypar0[2];

       }

# Set the initial value of count array.

      for (i=1;i<=ntotal;i++)

       {

       pair00[i]=0;

       }

# The parameter for the first read content in cell.

     nread=0;

     nstep=0;

     dist=deltar*ntotal;

print "Now, we will calculate the pair correlation function between ",atoma," and ", atomb,".";

print "The cell size is x=",cellx,", y=",celly,"and z=",cellz,".";

print "The dynamic steps are from ",stepmin,"to",stepmax,".";

print "The inteval of distance (delta r, in Angstrom) is,"deltar",", "and the total step is ",ntotal,".";

print "If the above parameters have some thing wrong, please change your file 'pait.txt'to correct them."

     }

/STEP/{

# Determine the atoms number from the first STEP.

     if ( nread == 0 )

        {

        natoma=0;

        natomb=0;

        for (i=1;i<=100000;i++)

          {

          getline;

          if ( NF == 1 ) { nread = $1; break;}

          if ( $1 == atoma ) natoma=natoma+1;

          if ( $1 == atomb ) natomb=natomb+1;

          }

        }

# The main program to count the distance.

     if ( ($2 >= stepmin) && ($2 <= stepmax) )

      {

     na=0;

     nb=0;

     nstep=nstep+1;

     printf "%10s",nstep;

     for (i=1;i<=nread;i++)

      {

     getline;

       if ( $1 == atoma )

          {

          na=na+1;

          ax[na]=$2;

          ay[na]=$3;

          az[na]=$4;

          }

       if ( $1 == atomb )

          {

          nb=nb+1;

          bx[nb]=$2;

          by[nb]=$3;

          bz[nb]=$4;

          }

         }

   for (ia=1;ia<=na;ia++)

       {

   for (ix=-1;ix<=1;ix++)

       {

   for (iy=-1;iy<=1;iy++)

       {

   for (iz=-1;iz<=1;iz++)

       {

   for (ib=1;ib<=nb;ib++)

       {

        Rtest=int(sqrt((ax[ia]-bx[ib]-cellx*ix)^2+(ay[ia]-by[ib]-celly*iy)^2+(az[ia]-bz[ib]-cellz*iz)^2)/deltar+0.5) ;

        if ( Rtest <= ntotal ) pair00[Rtest]=pair00[Rtest]+1;

       }

       }

       }

       }

      }

    }}

END{

    ncoord=0;

    Rho0=atomb/cellx/celly/cellz;

print "The below is the pair correlation function between ",atoma," and ", atomb,"." > "Pair_"atoma"_"atomb".txt";

print "The cell size is x=",cellx,", y=",celly,"and z=",cellz,"." > "Pair_"atoma"_"atomb".txt";

print "The dynamic steps are from ",stepmin,"to",stepmax,"." > "Pair_"atoma"_"atomb".txt";

print "The inteval of distance (delta r, in Angstrom) is,"deltar",", "and the total step is ",ntotal,"." > "Pair_"atoma"_"atomb".txt";

print "If the above parameters have some thing wrong, please change your file 'pait.txt'to correct them." > "Pair_"atoma"_"atomb".txt";

    for (i=1;i<=ntotal;i++)

      {

    ncoord=ncoord+pair00[i];

    printf "%8.3f%10.3f%10.3f\n",i*deltar, \

        pair00[i]/(natoma*natomb/cellx/celly/cellz)/(4*3.1415926*(i*deltar)^2*deltar)/nstep,ncoord/natoma/nstep > "Pair_"atoma"_"atomb".txt";

      }

   }

import ij.*;

import ij.plugin.filter.PlugInFilter;

import ij.process.*;

import ij.gui.*;

import ij.measure.Calibration;

import java.awt.*;

import java.util.*;

public class Radial_Profile implements PlugInFilter {

        ImagePlus imp;

        boolean canceled=false;

        double X0;

        double Y0;

        double mR;

        Rectangle rct;

        int nBins=100;

        //static boolean doNormalize = true;

        static boolean useCalibration = false;

        public int setup(String arg, ImagePlus imp) {

                this.imp = imp;

                return DOES_ALL+NO_UNDO+ROI_REQUIRED;

        }

        public void run(ImageProcessor ip) {

                setXYcenter();

                IJ.makeOval((int)(X0-mR), (int)(Y0-mR), (int)(2*mR), (int)(2*mR));

                doDialog();

                IJ.makeOval((int)(X0-mR), (int)(Y0-mR), (int)(2*mR), (int)(2*mR));

                imp.startTiming();

                if (canceled) return;

                doRadialDistribution(ip);

        }

        private void setXYcenter() {

                rct = imp.getRoi().getBoundingRect();

                X0 = (double)rct.x+(double)rct.width/2;

                Y0 =  (double)rct.y+(double)rct.height/2;

                mR =  (rct.width+rct.height)/4.0;

        }

        private void doRadialDistribution(ImageProcessor ip) {

                nBins = (int) (3*mR/4);

                int thisBin;

                float[][] Accumulator = new float[2][nBins];

                double R;

                double xmin=X0-mR, xmax=X0+mR, ymin=Y0-mR, ymax=Y0+mR;

                for (double i=xmin; i
                        for (double j=ymin; j
                                R = Math.sqrt((i-X0)*(i-X0)+(j-Y0)*(j-Y0));

                                thisBin = (int) Math.floor((R/mR)*(double)nBins);

                                if (thisBin==0) thisBin=1;

                                thisBin=thisBin-1;

                                if (thisBin>nBins-1) thisBin=nBins-1;

                                Accumulator[0][thisBin]=Accumulator[0][thisBin]+1;

                                Accumulator[1][thisBin]=Accumulator[1][thisBin]+ip.getPixelValue((int)i,(int)j);

                        }

                }

                Calibration cal = imp.getCalibration();

                if (cal.getUnit() == "pixel") useCalibration=false;

                Plot plot = null;

                if (useCalibration) {

                        for (int i=0; i
                                Accumulator[1][i] =  Accumulator[1][i] / Accumulator[0][i];

                                Accumulator[0][i] = (float)(cal.pixelWidth*mR*((double)(i+1)/nBins));

                        }

                        plot = new Plot("Radial Profile Plot", "Radius ["+cal.getUnits()+"]", "Normalized Integrated Intensity",  Accumulator[0], Accumulator[1]);

                } else {

                        for (int i=0; i
                                Accumulator[1][i] = Accumulator[1][i] / Accumulator[0][i];

                                Accumulator[0][i] = (float)(mR*((double)(i+1)/nBins));

                        }

                        plot = new Plot("Radial Profile Plot", "Radius [pixels]", "Normalized Integrated Intensity",  Accumulator[0], Accumulator[1]);

                }

                plot.show();

        }

        private void doDialog() {

                canceled=false;

                GenericDialog gd = new GenericDialog("Radial Distribution...", IJ.getInstance());

                gd.addNumericField("X center (pixels):",X0,2);

                gd.addNumericField("Y center (pixels):", Y0,2);

                gd.addNumericField("Radius (pixels):", mR,2);

                //gd.addCheckbox("Normalize", doNormalize);

                gd.addCheckbox("Use Spatial Calibration", useCalibration);

                gd.showDialog();

                if (gd.wasCanceled()) {

                        canceled = true;

                        return;

                }

                X0=gd.getNextNumber();

                Y0=gd.getNextNumber();

                mR=gd.getNextNumber();

                //doNormalize = gd.getNextBoolean();

                useCalibration = gd.getNextBoolean();

                if(gd.invalidNumber()) {

                        IJ.showMessage("Error", "Invalid input Number");

                        canceled=true;

                        return;

                }

        }

}

*******************************************************************
python(出处忘记了):
*******************************************************************

CODE:

#!/usr/bin/env python

# -*- coding: utf-8 -*-

# Generate RDF from PDB

# Copyright (C) 2011 Stas Bevc [email]stas.bevc@cmm.ki.si[/email]

# This program is free software: you can redistribute it and/or modify

# it under the terms of the GNU General Public License as published by

# the Free Software Foundation, either version 3 of the License, or

# (at your option) any later version.

# This program is distributed in the hope that it will be useful,

# but WITHOUT ANY WARRANTY; without even the implied warranty of

# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License

# along with this program.  If not, see <[url]http://www.gnu.org/licenses/[/url]>

from os import listdir

from math import sqrt, ceil, floor, pow, pi

# generate radial distribution function from

# particle positions stored in PDB files

def makeRDF(pdbdir, side, numbins, cm=False, numAT=0):

    maxbin = numbins # number of bins

    sideh = side/2.0

    dr = float(sideh/maxbin) # bin width

    hist = [0]*(maxbin+1) # initialize list

    rdf = {}

    count = 0

    pdbs = listdir(pdbdir) # list of files

    nstep = len(pdbs)

    print "Directory "+pdbdir+" has "+str(len(pdbs))+" files."

    # loop over pdb files

    for pdb in pdbs:

        if not pdb.endswith(".pdb"):

            nstep -= 1

            continue # skip other files

        count += 1

        print "Reading file ... "+pdb+" ("+ str(count) +"/"+ str(len(pdbs)) +")"

        # read atom coordinates from PDB

        atoms = []

        cmm = [0,0,0] # center of mass of molecule

        atc = 1 # atom counter

        lines = open(pdbdir+"/"+pdb)

        for line in lines:

            if line[0:4] != "ATOM":

                continue # skip other lines

            coords = map(float, (line[31:54]).split())

            if cm == True: # calculate center of mass

                cmm[0] += coords[0]

                cmm[1] += coords[1]

                cmm[2] += coords[2]

                if atc < numAT:

                    atc += 1

                else:

                    atc = 1

                    cmm[0] /= numAT

                    cmm[1] /= numAT

                    cmm[2] /= numAT

                    # fold coordinates

                    for i in range(3):

                        tmp = floor(cmm[i] * (1.0/side))

                        cmm[i] -= tmp * side

                    atoms.append((cmm[0],cmm[1],cmm[2]))

                    cmm = [0,0,0]

            else: # no cm calculation

                atoms.append((coords[0], coords[1], coords[2]))

        # loop over particle pairs

        npart = len(atoms)

        print " looping over particle pairs (" +str(npart)+ "^2) ... "

        for i in range(npart):

            xi = (atoms[i])[0]

            yi = (atoms[i])[1]

            zi = (atoms[i])[2]

            for j in range(i+1, npart):

                xx = xi - (atoms[j])[0]

                yy = yi - (atoms[j])[1]

                zz = zi - (atoms[j])[2]

                # minimum image

                if (xx < -sideh):   xx = xx + side

                if (xx > sideh):    xx = xx - side

                if (yy < -sideh):   yy = yy + side

                if (yy > sideh):    yy = yy - side

                if (zz < -sideh):   zz = zz + side

                if (zz > sideh):    zz = zz - side

                # distance between i and j

                rd  = xx * xx + yy * yy + zz * zz

                rij = sqrt(rd)

                bin = int(ceil(rij/dr)) # determine in which bin the distance falls

                if (bin <= maxbin):

                    hist[bin] += 1

    # normalize

    print "Normalizing ... "

    phi = npart/pow(side, 3.0) # number density (N*V)

    norm = 2.0 * pi * dr * phi * nstep * npart

    for i in range(1, maxbin+1):

        rrr = (i - 0.5) * dr

        val = hist[i]/ norm / ((rrr * rrr) + (dr * dr) / 12.0)

        rdf.update({rrr:val})

    return rdf

#-------------------------------------------------------------------#

# write RDF into file

boxsize = 36.845

numbins = 384 # number of bins

cm = True # calculate RDF from center of mass of molecule

numAT = 4 # number of atoms in molecule

pdbsdir = "./pdbs-ad40k-cg-ex-paral/" # directory with PDB files

outfile = "./rdf-ad40k-cg-ex-paral.out"

rdf = makeRDF(pdbsdir, boxsize, numbins, cm, numAT)

print "Writing output file ... " +outfile

outfile = open(outfile, "w")

for r in sorted(rdf.iterkeys()): # sort before writing into file

    outfile.write("%15.8g %15.8g\n"%(r, rdf[r]))

outfile.close()

******************************************************************
Fortran(参见ref.5)：
******************************************************************

CODE:

SUBROUTINE GR(NSWITCH)

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)

      PARAMETER(NM=40000,PI=3.141592653589793D0,NHIS=100)

      COMMON/LCS/X0(3,-2:2*NM),X(3,-2:2*NM,5),XIN(3,-2:2*NM),

     $XX0(3,-2:2*NM),XX(3,-2:2*NM,5),XXIN(3,-2:2*NM)

      COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM,NC,NN,MC

      COMMON/WALLS/HI(3,3),G(3,3),DH,AREA,VOLUME,SCM(3)

      COMMON/PBCS/HALF,PBCX,PBCY,PBCZ

        COMMON/GR_VAR/ NGR

        DIMENSION H(3,3),GG(0:NHIS),R(0:NHIS)

      EQUIVALENCE(X0(1,-2),H(1,1))

C   *****************************************************************

C      如何确定分子数密度：DEN_IDEAL 

C      取分子总数作为模拟盒中的数密度，可保证采样分子总数=总分子数?

C====================================================================

C        　N1=MOLSP+1

C      N2=MOLSP+NC

      DEN_IDEAL=MOLSP 

        G11=G(1,1)

      G22=G(2,2)

      G33=G(3,3)

      G12D=G(1,2)+G(2,1)

      G13D=G(1,3)+G(3,1)

      G23D=G(2,3)+G(3,2)

      IF(NSWITCH.EQ.0)THEN

          NGR=0

          DELR=HALF/NHIS

          DO I=1,NHIS

           GG(I)=0.D0 

           R(I)=0.D0 

          ENDDO 

      ELSE IF(NSWITCH.EQ.1)THEN

         NGR=NGR+1

       DO I=1,MOLSP-1

         DO J=I+1,MOLSP

C====================================================================

C     USE PBC IN X DIRECTION:  SUITABLE FOR PBCX=1

C                              NOT GREAT PROBLEM FOR PBCX=0 

C                              (THIS TIME USUALLY |DELTA X| < HALF)

C====================================================================

          XIJ=X0(1,I)-X0(1,J)

        IF(XIJ.GT.+HALF)XIJ=XIJ-PBCX

        IF(XIJ.LT.-HALF)XIJ=XIJ+PBCX

        YIJ=X0(2,I)-X0(2,J)

        IF(YIJ.GT.+HALF)YIJ=YIJ-PBCY

        IF(YIJ.LT.-HALF)YIJ=YIJ+PBCY

        ZIJ=X0(3,I)-X0(3,J)

        IF(ZIJ.GT.+HALF)ZIJ=ZIJ-PBCZ

        IF(ZIJ.LT.-HALF)ZIJ=ZIJ+PBCZ

        RSQ=XIJ*(G11*XIJ+G12D*YIJ+G13D*ZIJ)+

     $      YIJ*(G22*YIJ+G23D*ZIJ)+G33*ZIJ*ZIJ

          RRR=SQRT(RSQ)

          RRR=RRR/H(1,1)

C====================================================================

C      以上用数组G和H的结果与下同

C      RRR=SQRT(XIJ**2+YIJ**2+ZIJ**2)

C      G11=H(1,1)**2

C====================================================================

          IF(RRR.LT.HALF)THEN

           IG=INT(RRR/DELR)

           GG(IG)=GG(IG)+2

          ENDIF

       ENDDO

         ENDDO

      ELSE IF(NSWITCH.EQ.2)THEN

        DO I=1,NHIS

           R(I)=DELR*(I+0.5D0)

        ENDDO

        DO I=1,NHIS

           VB=(4.D0/3.D0)*PI*(((I+1)**3-I**3)*(DELR**3))

           GNID=VB*DEN_IDEAL

           GG(I)=GG(I)/(NGR*MOLSP*GNID)

        ENDDO

        OPEN(UNIT=31,FILE="GR.DAT")

        DO I=1,NHIS

           WRITE(31,*)R(I),GG(I)

        ENDDO

        CLOSE(31)

        ENDIF

        RETURN 

        END

**************************************************************************
awk(针对CPMD计算得出的TRAJEC.xyz)
**************************************************************************

CODE:

BEGIN{

#Set a pair of atoms you want to consider.

     atoma="O";

     atomb="O";

#Set the size of cell: x, y, z.

     cellx=11.3061;

     celly=11.3061;

     cellz=11.3061;

#Set the minimum and maxinum dynamics steps you want to calculated.

     stepmin=50000;

     stepmax=60000;

#Set the interval (delta r, in Angstrom) and the total steps of radius (r).

     deltar=0.05;

     ntotal=200;

# Read user's initial parameter from the file "pair.txt" in current directory.

     for (i=1;i<=100;i++)

       {

       getline mypar < "pair.txt";

       split(mypar,mypar0);

       if ( mypar0[1] == "atoma" ) atoma=mypar0[2];

       if ( mypar0[1] == "atomb" ) atomb=mypar0[2];

       if ( mypar0[1] == "cellx" ) cellx=mypar0[2];

       if ( mypar0[1] == "celly" ) celly=mypar0[2];

       if ( mypar0[1] == "cellz" ) cellz=mypar0[2];

       if ( mypar0[1] == "stepmin" ) stepmin=mypar0[2];

       if ( mypar0[1] == "stepmax" ) stepmax=mypar0[2];

       if ( mypar0[1] == "deltar" ) deltar=mypar0[2];

       if ( mypar0[1] == "ntotal" ) ntotal=mypar0[2];

       }

# Set the initial value of count array.

      for (i=1;i<=ntotal;i++)

       {

       pair00[i]=0;

       }

# The parameter for the first read content in cell. 

     nread=0;

     nstep=0;

     dist=deltar*ntotal;

print "Now, we will calculate the pair correlation function between ",atoma," and ", atomb,".";

print "The cell size is x=",cellx,", y=",celly,"and z=",cellz,".";

print "The dynamic steps are from ",stepmin,"to",stepmax,".";

print "The inteval of distance (delta r, in Angstrom) is,"deltar",", "and the total step is ",ntotal,".";

print "If the above parameters have some thing wrong, please change your file 'pait.txt'to correct them." 

     }

/STEP/{

# Determine the atoms number from the first STEP.

     if ( nread == 0 )

        {

        natoma=0;

        natomb=0;

        for (i=1;i<=100000;i++)

          {

          getline;

          if ( NF == 1 ) { nread = $1; break;}

          if ( $1 == atoma ) natoma=natoma+1;

          if ( $1 == atomb ) natomb=natomb+1;

          }

        }

# The main program to count the distance.

     if ( ($2 >= stepmin) && ($2 <= stepmax) )

      {

     na=0;

     nb=0;

     nstep=nstep+1;

     printf "%10s",nstep;

     for (i=1;i<=nread;i++)

      {

     getline;

       if ( $1 == atoma )

          {

          na=na+1;

          ax[na]=$2;

          ay[na]=$3;

          az[na]=$4;

          }

       if ( $1 == atomb )

          {

          nb=nb+1;

          bx[nb]=$2;

          by[nb]=$3;

          bz[nb]=$4;

          }

         }

   for (ia=1;ia<=na;ia++)

       {

   for (ix=-1;ix<=1;ix++)

       {

   for (iy=-1;iy<=1;iy++)

       {

   for (iz=-1;iz<=1;iz++)

       {

   for (ib=1;ib<=nb;ib++)

       {

        Rtest=int(sqrt((ax[ia]-bx[ib]-cellx*ix)^2+(ay[ia]-by[ib]-celly*iy)^2+(az[ia]-bz[ib]-cellz*iz)^2)/deltar+0.5) ;

        if ( Rtest <= ntotal ) pair00[Rtest]=pair00[Rtest]+1;

       }

       }

       }

       }

      }

    }}

END{

    ncoord=0;

    Rho0=atomb/cellx/celly/cellz;

print "The below is the pair correlation function between ",atoma," and ", atomb,"." > "Pair_"atoma"_"atomb".txt";

print "The cell size is x=",cellx,", y=",celly,"and z=",cellz,"." > "Pair_"atoma"_"atomb".txt";

print "The dynamic steps are from ",stepmin,"to",stepmax,"." > "Pair_"atoma"_"atomb".txt";

print

 "The inteval of distance (delta r, in Angstrom) is,"deltar",", "and the

 total step is ",ntotal,"." > "Pair_"atoma"_"atomb".txt";

print 

"If the above parameters have some thing wrong, please change your file 

'pait.txt'to correct them." > "Pair_"atoma"_"atomb".txt"; 

    for (i=1;i<=ntotal;i++)

      {

    ncoord=ncoord+pair00[i];

    printf "%8.3f%10.3f%10.3f\n",i*deltar, \

pair00[i]/(natoma*natomb/cellx/celly/cellz)/(4*3.1415926*(i*deltar)^2*deltar)/nstep,ncoord/natoma/nstep

 > "Pair_"atoma"_"atomb".txt";

      }

   }