[gmx-developers] reading velocity in gmx_msd.c
Mark Abraham
mark.abraham at anu.edu.au
Sat May 22 16:44:16 CEST 2010
----- Original Message -----
From: David van der Spoel <spoel at xray.bmc.uu.se>
Date: Sunday, May 23, 2010 0:08
Subject: Re: [gmx-developers] reading velocity in gmx_msd.c
To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> On 2010-05-22 15.56, Gaurav Goel wrote:
> >I've not used the debugger. I located the error point by
> writing to
> >stdout from different blocks in the code.
That works for a while, but what you really want to know is whether there's memory allocated for v, which a debugger will tell you right away. If you look back at the code for xa, there's an if statement regarding bMol that allocates xa or not. You probably didn't do that for va.
Mark
> did you allocate memory for the velocities?
> >
> >I'm using v. 4.0.5.
> >
> >Goal: Calculate viscosity with an analogous mean-squared displacement
> >method since it converges much faster. So, I need to practically
> >repeat the diffusivity calculation except that [r_i(t) -
> r_i(0)]**2 is
> >replaced by [r_ix(t)*p_iy(t) - r_ix(0)p_iy(0)]**2. This method
> >converges much faster... advantages and constraints on
> application of
> >this method are in Phys. Rev. A, 43(8), 4289 (1991).
> >
> >Diff against the code:
> >-----------------
> >1d0
> >< /*
> >75a75
> >> rvec **v0;
> >83c83
> >< typedef real t_calc_func(t_corr
> *,int,atom_id[],int,rvec[],rvec,bool,matrix);>---
> >>typedef real t_calc_func(t_corr
> *,int,atom_id[],int,rvec[],rvec[],rvec,bool,matrix);>111a112
> >> this-
> >v0 = NULL;
> >153c154
> >< for(i=0; (i<this->nrestart); i++)
> >---
> >> for(i=0; (i<this->nrestart); i++) {
> >154a156,157
> >> sfree(this->v0[i]);
> >> }
> >155a159
> >> sfree(this->v0);
> >195c199
> >< static void calc_corr(t_corr *this,int nr,int
> nx,atom_id index[],rvec xc[],
> >---
> >>static void calc_corr(t_corr *this,int nr,int nx,atom_id
> index[],rvec xc[],rvec vc[],
> >206a211
> >> memcpy(this->v0[this-
> >nlast],vc,this->ncoords*sizeof(vc[0]));
> >222c227
> >< g =
> calc1(this,nx,index,nx0,xc,dcom,bTen,mat);>---
> >> g =
> calc1(this,nx,index,nx0,xc,vc,dcom,bTen,mat);>235c240
> >< static real calc1_norm(t_corr *this,int nx,atom_id
> index[],int nx0,rvec xc[],
> >---
> >>static real calc1_norm(t_corr *this,int nx,atom_id index[],int
> nx0,rvec xc[],rvec vc[],
> >306c311
> >< static real calc_one_mw(t_corr *this,int ix,int
> nx0,rvec xc[],real *tm,
> >---
> >>static real calc_one_mw(t_corr *this,int ix,int nx0,rvec
> xc[],rvec vc[],real *tm,
> >350c355
> >< static real calc1_mw(t_corr *this,int nx,atom_id
> index[],int nx0,rvec xc[],
> >---
> >>static real calc1_mw(t_corr *this,int nx,atom_id index[],int
> nx0,rvec xc[],rvec vc[],
> >359c364
> >< g +=
> calc_one_mw(this,index[i],nx0,xc,&tm,dcom,bTen,mat);>---
> >> g +=
> calc_one_mw(this,index[i],nx0,xc,vc,&tm,dcom,bTen,mat);>429c434
> >< static real calc1_mol(t_corr *this,int nx,atom_id
> index[],int nx0,rvec xc[],
> >---
> >>static real calc1_mol(t_corr *this,int nx,atom_id index[],int
> nx0,rvec xc[],rvec vc[],
> >440c445
> >< g =
> calc_one_mw(this,i,nx0,xc,&mm,dcom,bTen,mat);>---
> >> g =
> calc_one_mw(this,i,nx0,xc,vc,&mm,dcom,bTen,mat);>529a535
> >>
> rvec *v[2],*va[2];
> >536c542,543
> >< natoms = read_first_x(&status,fn,&this-
> >t0,&(x[cur]),box);>---
> >> natoms = read_first_xv(&status,fn,&this-
> >t0,&(x[cur]),&(v[cur]),box);>>
> >545c552
> ><
> >---
> >> snew(v[prev],natoms);
> >549a557,558
> >> snew(va[0],this->ncoords);
> >> snew(va[1],this->ncoords);
> >553a563,564
> >> va[0] = v[0];
> >> va[1] = v[1];
> >575d585
> ><
> >577a588,591
> >>
> >> srenew(this->v0,this-
> >nrestart);>> snew(this-
> >v0[this->nrestart-1],this->ncoords);
> >>
> >615a630
> >>
> memcpy(va[prev],va[cur],this->ncoords*sizeof(va[prev][0]));
> >628c643,644
> ><
> >---
> >>
> >>
> >631a648
> >>
> >635c652,653
> ><
> calc_corr(this,i,gnx[i],index[i],xa[cur],bRmCOMM,com,calc1,bTen);>---
> >>
> >>
> calc_corr(this,i,gnx[i],index[i],xa[cur],va[cur],bRmCOMM,com,calc1,bTen);>641c659
> >< } while
> (read_next_x(status,&t,natoms,x[cur],box));>---
> >> } while
> (read_next_xv(status,&t,natoms,x[cur],v[cur],box));>905a924,925
> >>
> >>
> >946c966
> ><
> >---
> >>
> >--------------------------------------------------------
> >
> >-Gaurav
> >
> >On Fri, May 21, 2010 at 8:30 PM, Mark
> Abraham<mark.abraham at anu.edu.au> wrote:
> >>Did you use a debugger to find out why your segfault is arising?
> >>
> >>What GROMACS version are you using?
> >>
> >>A diff against the code with a few lines of context is much
> easier for us to understand, if we need to look at your code :-)
> >>
> >>I think I know what the problem is, but you haven't given
> enough information for me to be sure :-)
> >>
> >>Mark
> >>
> >>----- Original Message -----
> >>From: Gaurav Goel<gauravgoeluta at gmail.com>
> >>Date: Saturday, May 22, 2010 6:18
> >>Subject: Re: [gmx-developers] reading velocity in gmx_msd.c
> >>To: Discussion list for GROMACS development<gmx-
> developers at gromacs.org>>>
> >>>**************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.
> >>>>
> >>>--
> >>>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.
> >>--
> >>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.
> >>
>
>
> --
> 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.
More information about the gromacs.org_gmx-developers
mailing list