[gmx-users] Changing gmx_msd.c

[Mark Abraham]

Jennifer Williams wrote:
> Hi gmx_users,
> 
> I am trying to use NEMD to get transport diffusion coefficients.
> 
> I am adding an acceleration to my atoms in the z direction and then I 
> need to calculate this quantity:
> 
> 1/t <  ? [r(t) ? r(0)] >
> 
> My plan was to alter the code that calculates the MSD (gmx_msd.c located 
> in src/tools). This calculates;
> 
> 1/t < 1/N. ? [r(t) ? r(0)]^2 >
> 
> I presumed that all I needed to do was locate and remove the squared 
> function and the division by the number of molecules (1/N).
> 
> I would also change the factor for the conversion of nm^2/ps to 10e-5 
> cm^2/s since I would now be doing nm/ps to 10e-5 cm/s.

I'll have to guess you're correct, since I can't read the symbols you're 
intending.

> However my problem is that I have never worked with C++ so changing the 
> code involved some guess-work.

GROMACS is written in C, not C++, by the way :-)

> For the removal of the square I did;
> 
>     r2   += rv[m]*rv[m];   >                     r2   += rv[m];   (line 
> 252)
> 
> 
>         r2 += r*r;             >                     r2 += r;    (line 270)
> 
> to remove the division by the number of molecules (1/N), I removed line
> 
>   g/=nx;                 (line 279)
> 
> I realise that these changes are only valid for the first case CASE 
> (static_real_calc_norm), which I presumed was for the basic calculation 
> invoked by this command;

Assumption is a short route to madness when programming. You should look 
at the code from which the relevant function is called and see what 
preconditions are required for it. Stepping through the code with a 
debugger is probably a good way to understand the flow.

> g_msd ?f traj.xtc ?s md.tpr
> 
> i.e. I don?t plan on using ?pbc ?mol etc so I didn?t change the code for 
> these cases. (I also realise that unless I change the code, the 
> calculated error wont be valid-but I was going to worry about that later!)
> 
> However, the changes I made do not change the output from that of the 
> original MSD code. I am using the correct, changed executable and not an 
> old one (i.e I copied the compiled executable to gromacs/bin). So it 
> looks like I am not changing the correct part of the code, or perhaps 
> the command line g_msd ?f traj.xtc ?s md.tpr is not changing the case I 
> altered?

You can just run $builddir/src/tools/g_msd rather than "installing" a 
potentially problematic version. Adding a printf statement which is 
unique to your version will help you feel confident you're using the 
right version.

> Has anyone else changed the code for this purpose or would anyone be 
> able to give me some pointers on what I need to alter in the code. Today 
> is the first time I have ever looked at any C++ code (though I know some 
> fortran) so please be as specific as possible :)

There are easier things in life than modifying code written by someone 
else in an unfamiliar programming language (e.g. juggling knives while 
dictating a symphony) but you seem to be in the ballpark.

Mark