[gmx-developers] reading velocity in gmx_msd.c

[Gaurav Goel]

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]));
memcpy(xa[prev],xa[cur],this->ncoords*sizeof(xa[prev][0]));
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.
>