Re: [gmx-users] Changing gmx_msd.c

2009-07-18 Thread 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
___
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

Can't post? Read http://www.gromacs.org/mailing_lists/users.php


[gmx-users] Changing gmx_msd.c

2009-07-17 Thread Jennifer Williams

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.


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


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;


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?


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 :)



Much appreciated



--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.


___
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php