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;
                }
        }
}

#!/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()

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      如何確定分子數(shù)密度：DEN_IDEAL
C      取分子總數(shù)作為模擬盒中的數(shù)密度，可保證采樣分子總數(shù)=總分子數(shù)?
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      以上用數(shù)組G和H的結(jié)果與下同
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

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";
      }
   }