[gmx-developers] reading velocity in gmx_msd.c
Gaurav Goel
gauravgoeluta at gmail.com
Sat May 22 15:56:24 CEST 2010
I've not used the debugger. I located the error point by writing to
stdout from different blocks in the code.
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.
>
More information about the gromacs.org_gmx-developers
mailing list