[gmx-developers] g_anaeig generates wrong extreme projections when more than one eigenvector is analyzed

[Berk Hess]

Hi,

Next time please file a bugzilla.
I filed one now:
http://redmine.gromacs.org/issues/870
I posted a fix as well (there were more things to fix).
Please try it and report back on redmine.

Cheers,

Berk

On 01/25/2012 06:40 PM, Ivan Gushchin wrote:
> Hi,
>
> There seems to be a bug in g_anaeig (4.5.3) related to the 
> determination of the extreme (structural) projections of the 
> trajectory on the eigenvector, when more than one eigenvector is 
> analyzed. When I use different numbers of the eigenvectors for 
> analysis of the extreme projections, the generated fluctuations for 
> the same eigenvector differ in amplitude. The problem results, as I 
> understand, from the fact that the maximum and minimum projections 
> used for the generation of the output files (gmx_anaeig.c lines 
> 663-665) come from the last analyzed eigenvector (defined by -last).
>
> lines 624-636, determination of the extreme projections:
>
> for(v=0; v<noutvec_extr; v++) {
>     for(i=0; i<nframes; i++) {
>       if (inprod[v][i]<inprod[v][imin[v]])
>         imin[v]=i;
>       if (inprod[v][i]>inprod[v][imax[v]])
>         imax[v]=i;
>     }
> *min=inprod[v][imin[v]];
>     max=inprod[v][imax[v]];
> *fprintf(stderr,"%7d     %10.6f %10.1f %10.6f %10.1f\n",
>         eignr[outvec[v]]+1,
>         min,inprod[noutvec][imin[v]],max,inprod[noutvec][imax[v]]);
> }
>
> lines 663-665, generation of the output coordinates:
>
> xread[index[i]][d] =
>           (xav[i][d] + (*min**(nextr-frame-1)+*max**frame)/(nextr-1)
>           *eigvec[outvec[v]][i][d]/sqrtm[i]);
>
> As it can be seen, for all the eigenvectors the values of the extreme 
> projections used are the same, which is incorrect, as the vectors are 
> normalized but fluctuations along different vectors do differ in 
> amplitude. It can be healed by the following modification:
>
> xread[index[i]][d] =
>           (xav[i][d] + 
> (*inprod[v][imin[v]]**(nextr-frame-1)+*inprod[v][imax[v]]**frame)/(nextr-1)
>           *eigvec[outvec[v]][i][d]/sqrtm[i]);
>
> Could someone check if this is correct?
>
> Please excuse me if I post to the wrong mailing list.
>
> All the best,
> Ivan.
>
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20120130/45b09a6e/attachment.html>