[gmx-developers] reading velocity in gmx_msd.c

[Gaurav Goel]

**************Please ignore the previous mail... I accidentally
pressed enter before completing the mail***************


To be able to read velocities in g_msd I decided to do the following:

1. Add the following lines to src/gmxlib/trxio.c

int read_first_xv(int *status,char *fn,
                 real *t,rvec **x, rvec **v,matrix box)
{
 t_trxframe fr;

 read_first_frame(status,fn,&fr,TRX_NEED_X | TRX_NEED_V);
 if (*status >= nxframe) {
   nxframe = *status+1;
   srenew(xframe,nxframe);
 }
 xframe[*status] = fr;
 *t = xframe[*status].time;
 *x = xframe[*status].x;
 *v = xframe[*status].v;
 copy_mat(xframe[*status].box,box);

 return xframe[*status].natoms;
}


bool read_next_xv(int status,real *t, int natoms, rvec x[], rvec v[], matrix box
)
{
 bool bRet;

 xframe[status].x = x;
 xframe[status].v = v;
 bRet = read_next_frame(status,&xframe[status]);
 *t = xframe[status].time;
 copy_mat(xframe[status].box,box);

 return bRet;
}

2. Add the following lines to /include/statutil.h
extern int read_first_xv(int *status,char *fn,
                       real *t,rvec **x, rvec **v,matrix box);
extern bool read_next_xv(int status,real *t,int natoms,rvec x[], rvec v[],matrix
 box);

3. Make following changes to the gmx_msd.c subroutine:
/*define the followign pointers*/
rvec         *v[2],*va[2];

/*read x and v with these commands*/
natoms = read_first_xv(&status,fn,&this->t0,&(x[cur]),&(v[cur]),box);
while (read_next_xv(status,&t,natoms,x[cur],v[cur],box));

/*repeat the following commands as used for x*/
snew(v[prev],natoms);
  va[0] = v[0];
   va[1] = v[1];
srenew(this->v0,this->nrestart);
     snew(this->v0[this->nrestart-1],this->ncoords);

   if (bFirst) {

     memcpy(xa[prev],xa[cur],this->ncoords*sizeof(xa[prev][0]));
/*error point*/
memcpy(va[prev],va[cur],this->ncoords*sizeof(va[prev][0]));
bFirst = False;
}

Everything compiled all right. However, when I tried to use this new
utility I get a
'segmentation fault' error at
***memcpy(va[prev],va[cur],this->ncoords*sizeof(va[prev][0]));****


Am I not reading the velocities properly? Is my strategy of repeating
everything that was done for reading positions for velocities (except
for doing things like remove pbc) not correct?

I will greatly appreciate your help.

Thanks,
Gaurav

On Tue, May 18, 2010 at 5:16 PM, David van der Spoel
<spoel at xray.bmc.uu.se> wrote:
> On 2010-05-18 23.10, Gaurav Goel wrote:
>>
>> I am modifying the 'g_msd' utility to be able to calculate viscosities
>> using the corresponding Einstein relationship according to equation
>> 3.14 in Helfand E., Physical Review, 119, 1, 1960. The equation is
>> similar to the mean-square displacement equation for self-diffusion
>> coefficient with the position vector 'r_i ' replaced by 'r_ix p_iy'
>> where p_iy is the particle momenta in y-direction.
>>
>> So basically I need to be able to read both position and velocities
>> from the trajectory. Presently, gmx_msd.c uses 'read_first_x' and
>> 'read_next_x' to read the coordinates. Can you please comment on how I
>> can modify this code to read the corresponding velocities?
>>
>> Any help is greatly appreciated.
>>
>> thanks,
>> Gaurav
>
> From include/statutil.h:
>
> extern bool read_next_frame(const output_env_t oenv,int status,t_trxframe
> *fr);
>  /* Reads the next frame which is in accordance with fr->flags.
>   * Returns TRUE when succeeded, FALSE otherwise.
>   */
>
> Do also checkout
> g_tcaf -h
>
>
> --
> David van der Spoel, Ph.D., Professor of Biology
> Dept. of Cell & Molec. Biol., Uppsala University.
> Box 596, 75124 Uppsala, Sweden. Phone:  +46184714205.
> spoel at xray.bmc.uu.se    http://folding.bmc.uu.se
> --
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www interface
> or send it to gmx-developers-request at gromacs.org.
>