[CPMD-list] Pair correlation function calculation

[Axel Kohlmeyer]

On Sat, 30 Apr 2005, Younes Ansari wrote:

YA> Dear CPMD-list:

YA>   I have written a pair correlation function program using CPMD trajec.xyz by MATLAB6P:

YA> 1- I told it to input the No. of (1-iterations 2-timestep 3-volume
YA> 4-the distance you want to calculate from 0 to the distance 5-the
YA> Delta:  so

YA> if r =1 the program calculates r+delta to r-deta and if there is a R
YA> between this distance it will add LL=LL+1 else LL remains 0. )

YA> 2- I have calculated the distances for each iteration=R11 to Rij

YA> 3- calculations of RHO: (V/Z)*6.02E-7

YA> 4- my g(r)=(4/3)*3.141595*RHO*((1-delta)^3))*(LL)/(average)


dear younes ansari,

your scheme seems a little complicated and - as far as i can tell - 
does not consider the periodic boundary conditions. there are many
ways how you can do this. here is what i prefer to do:

lets assume you want to calculate the oxygen-hydrogen radial distribution
function for a bulk water system with 32 water molecules. lets also 
assume, that in your input file the 32 oxygen coordinates come first 
and then the 64 hydrogens, that means that the coordinates will be
in the same order in your trajectory file. and finally you need to 
know the dimensions of your supercell for the PBC treatment: let's 
assume a cubic box with 10.0 angstrom sidelength, but the following
scheme will work with any orthorhombic supercell.

now you you read in the trajectory step by step and for each step
you loop over the oxygen indices and then over the hydrogen indices
to collect the oxygen-hydrogen distances into a histogram with delta(r) 
0.1 angstrom:


do i=1,32
   do j=33,96
      call rdftohist(x(i),y(i),z(i),x(j),y(j),z(j),hist,0.1,boxx,boxy,boxz)
   end do
end do

in the rdftohist subroutine
you now add or subtract boxx so many times from dx=xi-xj that its value
is between -boxx/2 and boxx/2 and the same for dy and dz.
now you compute the minimum image distance from sqrt(dx*dx+dy*dy+dz*dz)
divide by delta add 0.5 and make it integer. this number is the histogram
index, that you have to increment by one to count this pair.

this you repeat for each time step and divide by the number of time steps.
each histogram entry i now represents the number of hydrogens in the 
distance from delta*(i-1) to delta*i. the final rdf is the ratio of this
quantity to a ideal gas distribution in the same volume. the normalization
coefficient is thus the volume of the reference particle divided by the 
volume of each sphere shell of the histogram

particle_vol = boxx * boxy * boxz / 32;
do i=1,nhist
   slice_vol = 4.0 / 3.0 * PI * ((delta*i)**3 - ((delta*(i-1))**3))
   g(i) = hist(i)*particle_vol/slice_vol
end do

to have g(oo) -> 1 you can additionally multiply by the ratio 
between the reference and pair particles (= 0.5 in this case).

YA> So I want to Know if You can send me a simple program of calculating
YA> g(r) using trajec.xyz. I can tell the program to read data from which


YA> columns and line of the file trajec.xyz. just I want a SCHEMATIC type
YA> of the g(r) calcultion. I have seen the books on the CPMD refrences
YA> [1,4] but there is not enough details to get it .  I want to know how
YA> to normalize my data.

i beg to differ: equation 2.94 in allen and tildesley describes _exactly_
the algorithm and the following section give a few hints on how to 
implement it. and section 6.2 outlines in detail the same algorithm as 
above for a monoatomic liquid (where you can save time by taking advantage 
of the fact, that the pair A-B is the same a B-A).

algorithm 7 in frenkel and smit is another version of what you want to 
calculate, also with example code.

regards,
	axel kohlmeyer.

YA> 
YA> 
YA> __________________________________________________
YA> Do You Yahoo!?
YA> Tired of spam?  Yahoo! Mail has the best spam protection around 
YA> http://mail.yahoo.com 

-- 

=======================================================================
Dr. Axel Kohlmeyer   e-mail: axel.kohlmeyer at theochem.ruhr-uni-bochum.de
Lehrstuhl fuer Theoretische Chemie          Phone: ++49 (0)234/32-26673
Ruhr-Universitaet Bochum - NC 03/53         Fax:   ++49 (0)234/32-14045
D-44780 Bochum  http://www.theochem.ruhr-uni-bochum.de/~axel.kohlmeyer/
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.