[gmx-users] number of solvent molecule

[chris.neale at utoronto.ca]

> thanks for ur reply. But could you please tell me the
> command to integrate the rdf to first maximum. I
> searched for it and did not find any thing and why I
> double the value?.

There is no command for integration. Mark's suggestion  
http://www.gromacs.org/pipermail/gmx-users/2006-November/024896.html
is the way to do it with g_rdf. I have never used g_hbond for this  
purpose, but Erik's suggestion  
http://www.gromacs.org/pipermail/gmx-users/2006-November/024904.html
sounds good for that.

If you do choose to use g_rdf, then use a spreadsheet to do the  
integration. Be aware though that what you want is not a simple  
cumulative value of g(r), but of N(dr) as I will outline below.

g(r) = N(dr)/V(dr) / N/V

Therefore to get the integral of N(dr) to a specified value, you will  
need to integrate g(r) * N/V * V(dr). You can easily get N from the  
number of waters in your system. V(dr) can be calculated in  
spreadsheet format by using a formula. In MS excel for example, you  
need many columns:

The first column is your r value and the second is your g(r), both  
columns 1 and 2 from the g_rdf output .xvg file. Columns 3 and 4  
should be defined as suggested below

r        g(r)             sphereVolume                     shellVolume
0.001    g(r)|dr=0.002    =4/3*3.1416*power(A1+0.001,3)    =C1
0.003    g(r)|dr=0.004    =4/3*3.1416*power(A2+0.001,3)    =C2-C1
0.005    g(r)|dr=0.006    =4/3*3.1416*power(A3+0.001,3)    =C3-C2

So your V(dr) value can come from the shellVolume column.
Copying and pasting the formulas in columns 3 and 4 will reproduce  
them as desired (although you need to create rows 1 and 2 and then  
copy row 2 on down). Also, if you don't use dr=0.002, then you will  
need to make a few changes in the formulas. Also note that r in .xvg  
is the central point along the r bin so that integrating to the trird  
row above would be out to 0.006nm not 0.005nm.

If you use the -cn option then the g(r) value that you should use is  
the one from the column at your desired r value. If not, then you need  
to sum over g(r) yourself.

The V (total volume) value is a little trickier. This entire procedure  
assumes that you can use a single value for V even though it changes  
every timestep. If your system is well behaved, then the difference  
should be tiny. I have previously posted a code snippet for g_vol to  
get the time averaged system volume  
http://www.gromacs.org/pipermail/gmx-users/2006-November/024694.html

If you end up using this menthod and also g_hbond, please post your  
general findings about their equivalence back to this list.