[gmx-users] number of solvent molecule
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Wed Nov 29 18:51:47 CET 2006
> 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.
More information about the gromacs.org_gmx-users
mailing list