Integration of g(r)

User discussions
1

GROMACS version: 2018
GROMACS modification: No
Here post your question

Hi everyone,

Recently, I am confused by the result of integration of g®. In the manul, the definition of g®_AB is g®_AB = pho®/ pho_local, the r_max for pho_local is the half of box.

Based on the definition, the integration from zero to rmax of g®_AB*pho_local will be the number of B (smaller than total number of B) in the sphere with radius r_max.

However, when I calculated the integration, I found that the result is bigger than the total number of B.

For example, if there are 26 A monomers and 26 B monomers in the system, the integration of g_ABpho_local will smaller than 26. But I got the value larger than 46.
The pho_local = average_number/(4.0/3pirmaxrmax*rmax), here average_number is accumulate number of B at distance rmax got from the output file using gmx rdf -cn

the integration script i used:

double dr = rr[1]-rr[0];
sumb = 0;
double psv = 0.0,rri;
for(i = 0; i < ri; i ++){
    rri = (i + 0.5) * dr;
    double sv = (4.0 / 3.0) * PI * rri * rri * rri;
    double dsv = sv-psv;
    sumb += (gr[i])*(dsv);
    psv = sv;
    cout<<"sumb = "<<sumb<<endl;
}
double pho_local = 24.8/((4.0 / 3.0) * PI * rri * rri * rri);
cout<<"total monomer number ~: "<<sumb*pho_local<<endl;;

I do not find anything wrong in my script. So could anyone help me to explain why the result is unexpected. Is there some misunderstanding of the definition of g® that GROMACS used to analyze data?

Thanks in advance

Best

isimuly

2

Hi isimuly,

purely following the math of (5.388) in the manual (http://manual.gromacs.org/2021-beta1/manual-2021-beta1.pdf) and using A=B you would obtain the number of particles from the distribution function by integrating the following way:

\langle \rho_{local}\rangle \int dr 4\pi r^2 g_{AA}(r) = 1/N_A \int d r \sum_{i\in A}\sum_{j\in A}\delta(r_{ij}-r) = N_A^2/N_A = N_A

(the integral on the right hand side over the delta function just counts the A,A pairs)

I agree that it’s a bit counter-intuitive to use 4\pi r^2 instead of the difference of the sphere volumes; it might help to think in terms of the continuous integral first and then approximate. There, you see that you obtain the factor by integrating in spherical coordinates over the complete angle range for \phi, \psi for a volume element r^2 \sin\phi \,\mathrm{d}\phi \,\mathrm{d}\psi \,\mathrm{d}r