[gmx-developers] reading velocity in gmx_msd.c
David van der Spoel
spoel at xray.bmc.uu.se
Sat May 22 16:08:02 CEST 2010
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.
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
More information about the gromacs.org_gmx-developers
mailing list