[gmx-developers] reading velocity in gmx_msd.c
David van der Spoel
spoel at xray.bmc.uu.se
Mon May 24 19:38:19 CEST 2010
On 5/24/10 5:07 PM, Gaurav Goel wrote:
> (gdb) where
> #0 0x0000003e2b482b6b in memcpy () from /lib64/libc.so.6
> #1 0x0000000000403992 in corr_loop ()
> #2 0x00000000004066df in do_corr ()
> #3 0x0000000000407f84 in gmx_visco ()
> #4 0x000000000040226b in main ()
>
>
You need to
cd src/tools
make clean
make CFLAGS=-g install
in order to be able to debug.
> -Gaurav
>
> On Mon, May 24, 2010 at 11:01 AM, David van der Spoel
> <spoel at xray.bmc.uu.se> wrote:
>
>> On 2010-05-24 16.47, Gaurav Goel wrote:
>>
>>> I had allocated memory for both v and va using
>>> snew(v[prev],natoms);
>>>
>>> and
>>> snew(va[0],this->ncoords);
>>> snew(va[1],this->ncoords);
>>>
>>> Results from the DEBUG script:
>>>
>>> $: more DEBUG
>>> #!/bin/bash -f
>>> xterm -e gdb -x gdb_cmds
>>> /home/gxg133/Download/gromacs-4.0.5/src/tools/g_visco
>>>
>>> $: more gdb_cmds
>>> break _gmx_error
>>> break gmx_visco
>>> run
>>>
>>> $: ./DEBUG
>>>
>>> Breakpoint 1 at 0x411670
>>> Breakpoint 2 at 0x4079f0
>>>
>>> Breakpoint 2, 0x00000000004079f0 in gmx_visco ()
>>> Missing separate debuginfos, use: debuginfo-install
>>> e2fsprogs-libs-1.41.4-12.fc11.x86_64 glibc-2.10.2-1.x86_64
>>> libICE-1.0.4-7.fc11.x86_64 libSM-1.1.0-4.fc11.x86_64
>>> libX11-1.2.2-1.fc11.x86_64 libXau-1.0.4-5.fc11.x86_64
>>> libxcb-1.2-5.fc11.x86_64
>>>
>>> (gdb) n
>>> Single stepping until exit from function gmx_visco,
>>> which has no line number information.
>>> :-) G R O M A C S (-:
>>>
>>> GRoups of Organic Molecules in ACtion for Science
>>>
>>> :-) VERSION 4.0.5 (-
>>> .....
>>> .....
>>>
>>> Group 0 ( System) has 2160 elements
>>> Group 1 ( Si1) has 720 elements
>>> Group 2 ( O1) has 1440 elements
>>> Select a group: 1
>>> Selected 1: 'Si1'
>>> Reading frame 0 time 0.000
>>> Program received signal SIGSEGV, Segmentation fault.
>>> 0x0000003e2b482b6b in memcpy () from /lib64/libc.so.6
>>> ----------------------
>>>
>>> So it seems the error is happening in "memcpy()". Any suggestions on
>>> how I proceed from here?
>>>
>>>
>> In an interactive debugger session you go up up up until you can see where
>> it is going wrong (or type "where").
>>
>>
>>
>>> -Gaurav
>>>
>>> On Sat, May 22, 2010 at 10:44 AM, Mark Abraham<mark.abraham at anu.edu.au>
>>> wrote:
>>>
>>>> ----- 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.
>>>>>
>>>> --
>>>> 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.
>>
>>
--
David.
________________________________________________________________________
David van der Spoel, PhD, Professor of Biology
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596, 75124 Uppsala, Sweden
phone: 46 18 471 4205 fax: 46 18 511 755
spoel at xray.bmc.uu.se spoel at gromacs.org http://xray.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
More information about the gromacs.org_gmx-developers
mailing list