[gmx-developers] reading velocity in gmx_msd.c
Gaurav Goel
gauravgoeluta at gmail.com
Mon May 24 16:47:07 CEST 2010
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?
-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.
>
More information about the gromacs.org_gmx-developers
mailing list