I need to calculate the radial distribution function for my system which includes 50 molecules. Here is my code line to generate rdf in gromacs:
gmx rdf -n Index.ndx -f mdrun1.xtc -s mdrun1.tpr -selrpos 'whole_res_com' -ref "AsphB" -sel "AsphB" -o rdf_com_20.xvg
I am wondering if I should choose all 50 molecules in both -ref and -sel sections?
Or should I consider molecule 1 as -ref and the rest as -sel and do that multiple times for all molecules and then get average overall? If yes, what type of average should that be?
Hi Ali! When you specify -ref and -sel in gromacs you specify all molecules which belong to the same group. If you have 1000 water molecules in your -ref then the software will give you average values for all of them (you do not need to do 1000 calculation for each molecule, this is rather a madness). The same applies for -sel. If you want to see rdfs for different molecules then you have to specify them in another calculation. I can add that older versions of gromacs could handle 1 reference versus 10 selected. So you could have several rdfs computed simultaneously. If you are seeking such an option then you have to install the older version too (like 4.6.7) and create a tpr-file for that version. Then you will have nice options for several groups in selection.
Hi Inna Ermilova, thank you so much for your reply.
I think I got your point. So if I have 1000 water molecule and I want to get the rdf for water molecules with themselves. I should prepare an index including all 1000 water molecules and consider that for both -ref and -sel. Is that what you meant?
Because I did that; however, there is an issue about that! I will have a huge pick at a very short distance that does not make sense! I think the reason for that is the consideration of each molecule with itself. For example, the distance between water molecule residue 1 in -ref with water molecule residue 1 in -sel! and considering that for all 1000 molecules will result in a huge pick in a very short distance. I tried to use -cut option but first I want to make sure if I am doing right!
I wonder if you have any idea. Your time and help are greatly appreciated.
Hi Ali! A huge peak at the short distance can mean that there is a hydrogen bonding between water molecules. You have to normalize your rdf. If you use the cut option you shall think the possible distances for hydrogen bonding. If you exclude those distances you will get the wrong distribution. And , yes, ref and sel should be the same if you compute for the same group. As I said, if you aim to get a normalized rdf which smoothly goes to 1 after 3.5 Å then you probably would like to use and older version of gromacs for it. I haven't seen -fade implemented in newer versions.
Hi Inna Ermilova, Thank you so much. It is very good information. I will definitely use your instruction for normalizing the RDF after I make sure I am not cutting off anything by mistake.
So, my molecule is not water actually and it does not have the potential to make a hydrogen bond. I have attached my molecule structure, dimension and the result of rdf if I consider all 50 molecules in both -ref and -sel. As you can see there is a huge pick before 0.2 nm (pulse a few more). Considering my molecule dimension and structure, I am wondering if those mentioned picks were generated mistakenly. I would appreciate you if you take look at the provided information and let me know what you think.
Please feel free to ask me for future information if you need. Your help means a lot to me. Thank you so much for your time.
HI Ali! I can not tell you if it is a mistake or not if I have no information about how you were doing this calculation, i.e. with what settings. For instance, if you didn't define more precise groups of atoms for calculating the rdfs, but just took groups for your molecules and used atom-atom rdf then this is a correct picture for the calculation. There are very many options which you can put in such an analysis. I would suggest you to read the manual for this routine and do the calculation more properly, because it doesn't look like you know what you defined there.
Hi, Inna Ermilova, thanks again for your reply. I would definitely reread the manual. But here is what I have done and know. My code line is:
gmx rdf -n Index.ndx -f mdrun1.xtc -s mdrun1.tpr -selrpos res_cog -ref AsphC -sel AsphC -o rdf_cog_totallast20.xvg -tu ns -dt 0.5 -b 100 -e 120
I intended to consider the centre of geometry of the whole molecule. my molecule has 112 atoms and I have 50 molecules. you can see the index that I used in the attached picture. I did it for the last 20 ns of simulation (the whole simulation time is 120ns) and I used the outputs every 500 ps (while my output store every 10 ps)
Hi Ali! So you calculate the rdfs between the centers of geometries of your residues. The output is fine for what you did. If you want to have a more clear result you shall make an index file which shall contain groups of atoms between which you want to compute your rdfs, then you shall select selrpos atom instead of res_cog. Or if you want to compute it for the center of mass of molecules then you shall write mol_com. You shall specify seltype here in the same fashion then, which you didn't do in that line.
Hi Inna Ermilova, Thank you so much for being so helpful and informative. I really appreciate you.
- Badges
- Science topic
- Similar topics
- Chemistry
- Theoretical Chemistry
- Computational Chemistry
- Molecular Dynamics