[gmx-developers] Coordinate scaling in pdbio.c

[Mark Abraham]

On 11/04/12, Erik Marklund <erikm at xray.bmc.uu.se> wrote:

> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 11 apr 2012 kl. 01.38 skrev Mark Abraham:
> 
> 
> 
> > 
> > On 11/04/2012 5:21 AM, David van der Spoel wrote:
> > 
> > 
> > > On 2012-04-10 21:10, Anton Feenstra wrote:
> > > 
> > 
> > 
> > > 
> > > 
> > > > On 09/04/12 21:47, Justin A. Lemkul wrote:
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > Radhakrishna Bettadapura wrote:
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > All,
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > I'm trying to figure out why a pdb input to pdb2gmx results in
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > coordinates that are scaled by a factor of 0.1. That is, if a line in
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > the PDB file reads
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > ATOM 1 N ALA B 3 -1.221 20.481 12.450 1.00 34.64 N
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > then the corresponding line in the .gro file reads
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 3ALA N 1 -0.122 2.048 1.245
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > The culprit seems to be the following lines in pdbio.c:
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > static int read_atom(t_symtab *symtab,
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > char line[],int type,int natom,
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > t_atoms *atoms,rvec x[],int chainnum,gmx_bool bChange)
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > {
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > // do stuff
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > x[natom][XX]=strtod(xc,NULL)*0.1; /*all coordinates scaled by 0.1...
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > why?*/
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > x[natom][YY]=strtod(yc,NULL)*0.1;
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > x[natom][ZZ]=strtod(zc,NULL)*0.1;
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > //do more stuff
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > }
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > There's also a line in the output function that multiplies all
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > incoming coordinates by 10 before writing it to output:
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > void write_pdbfile_indexed(FILE *out,const char *title,
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > t_atoms *atoms,rvec x[],
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > int ePBC,matrix box,char chainid,
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > int model_nr, atom_id nindex, atom_id index[],
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > gmx_conect conect, gmx_bool bTerSepChains)
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > {
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > //do stuff...
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > fprintf(out,pdbform,pdbtp[type],(i+1)%100000,nm,resnm,ch,resnr,
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > (resic == '\0') ? ' ' : resic,
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 10*x[i][XX],10*x[i][YY],10*x[i][ZZ],occup,bfac,atoms->atom[i].elem);
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > // do more stuff...
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > }
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > Can someone tell me why this coordinate scaling occurs? And why the
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > scale factor is a single hard-coded number than, say, something that
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > depends on the input?
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > PDB coordinates are in Angstrom, while .gro coordinates are in nm. Hence
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > > > the conversion factor of 0.1.
> > > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > What happened to A2NM and NM2A?
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > In my (4.0.5, yes that is very old) includes/physics.h I still see:
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > #define A2NM (ANGSTROM/NANO) /* NANO */
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > #define NM2A (NANO/ANGSTROM) /* 10.0 */
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > Why aren't these still used for that?
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > 
> > 
> > 
> > > 
> > > 
> > > > 
> > > > 
> > > 
> > 
> > 
> > > At some stage macros were removed from the code, even though these seem rather harmless. I agree that it is confusing to hard code these numbers. If we cannot use macro's like this we should probably replace them by
> > > 
> > 
> > 
> > > static const real ANGSTROM=1e-10;
> > > 
> > 
> > 
> > > static const real NANO=1e-9;
> > > 
> > 
> > 
> > > 
> > > 
> > 
> > 
> > > etc.
> > > 
> > 
> > 
> > > 
> > > 
> > 
> > 
> > > Comments?
> > > 
> > 
> > include/physics.h still has these macros. IMO, hard-coded constants are a greater evil than macros to prevent that, though I expect we will transition to const values at some stage soon.
> > 
> > 
> 
> 
> 
> But are these the kind of macros we want to avoid? Aren't function-like macros the ones to kill in the first place?
> 
> 
> 

 They're taking no parameters, so they're hardly function-like. The code fragments above are compiled into constants by the pre-processor. A macro that is used as
 
dist_in_nm = dist_in_angstrom * A2NM;
 
is much less evil than a macro 
 
dist_in_nm = A2NM(dist_in_angstrom);
 
and physics.h has only the first type. Inlined functions of the second type would be best of all, of course. 
 
Mark 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20120411/e3b61eff/attachment.html>