3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
52 #include "gmx_fatal.h"
58 #include "eigensolver.h"
62 static real
find_pdb_bfac(t_atoms
*atoms
,t_resinfo
*ri
,char *atomnm
)
67 strcpy(rresnm
,*ri
->name
);
69 for(i
=0; (i
<atoms
->nr
); i
++) {
70 if ((ri
->nr
== atoms
->resinfo
[atoms
->atom
[i
].resind
].nr
) &&
71 (ri
->ic
== atoms
->resinfo
[atoms
->atom
[i
].resind
].ic
) &&
72 (strcmp(*atoms
->resinfo
[atoms
->atom
[i
].resind
].name
,rresnm
) == 0) &&
73 (strstr(*atoms
->atomname
[i
],atomnm
) != NULL
))
77 fprintf(stderr
,"\rCan not find %s%d-%s in pdbfile\n",
78 rresnm
,ri
->nr
,atomnm
);
82 return atoms
->pdbinfo
[i
].bfac
;
85 void correlate_aniso(const char *fn
,t_atoms
*ref
,t_atoms
*calc
,
86 const output_env_t oenv
)
91 fp
= xvgropen(fn
,"Correlation between X-Ray and Computed Uij","X-Ray",
93 for(i
=0; (i
<ref
->nr
); i
++) {
94 if (ref
->pdbinfo
[i
].bAnisotropic
) {
95 for(j
=U11
; (j
<=U23
); j
++)
96 fprintf(fp
,"%10d %10d\n",ref
->pdbinfo
[i
].uij
[j
],calc
->pdbinfo
[i
].uij
[j
]);
102 static void average_residues(double f
[],double **U
,int uind
,
103 int isize
,atom_id index
[],real w_rls
[],
112 for(i
=0; i
<isize
; i
++) {
113 av
+= w_rls
[index
[i
]]*(f
!= NULL
? f
[i
] : U
[i
][uind
]);
114 m
+= w_rls
[index
[i
]];
116 atoms
->atom
[index
[i
]].resind
!=atoms
->atom
[index
[i
+1]].resind
) {
119 for(j
=start
; j
<=i
; j
++)
122 for(j
=start
; j
<=i
; j
++)
132 void print_dir(FILE *fp
,real
*Uaver
)
134 real eigvec
[DIM
*DIM
];
139 fprintf(fp
,"MSF X Y Z\n");
142 fprintf(fp
," %c ",'X'+d
-XX
);
144 fprintf(fp
," %9.2e",Uaver
[3*m
+d
]);
145 fprintf(fp
,"%s\n",m
==DIM
? " (nm^2)" : "");
148 for(m
=0; m
<DIM
*DIM
; m
++)
152 eigensolver(tmp
,DIM
,0,DIM
,eigval
,eigvec
);
154 fprintf(fp
,"\n Eigenvectors\n\n");
155 fprintf(fp
,"Eigv %-8.2e %-8.2e %-8.2e (nm^2)\n\n",
156 eigval
[2],eigval
[1],eigval
[0]);
159 fprintf(fp
," %c ",'X'+d
-XX
);
160 for(m
=DIM
-1; m
>=0; m
--)
161 fprintf(fp
,"%7.4f ",eigvec
[3*m
+d
]);
166 int gmx_rmsf(int argc
,char *argv
[])
168 const char *desc
[] = {
169 "g_rmsf computes the root mean square fluctuation (RMSF, i.e. standard ",
170 "deviation) of atomic positions ",
171 "after (optionally) fitting to a reference frame.[PAR]",
172 "With option [TT]-oq[tt] the RMSF values are converted to B-factor",
173 "values, which are written to a pdb file with the coordinates, of the",
174 "structure file, or of a pdb file when [TT]-q[tt] is specified.",
175 "Option [TT]-ox[tt] writes the B-factors to a file with the average",
177 "With the option [TT]-od[tt] the root mean square deviation with",
178 "respect to the reference structure is calculated.[PAR]",
179 "With the option [TT]aniso[tt] g_rmsf will compute anisotropic",
180 "temperature factors and then it will also output average coordinates",
181 "and a pdb file with ANISOU records (corresonding to the [TT]-oq[tt]",
182 "or [TT]-ox[tt] option). Please note that the U values",
183 "are orientation dependent, so before comparison with experimental data",
184 "you should verify that you fit to the experimental coordinates.[PAR]",
185 "When a pdb input file is passed to the program and the [TT]-aniso[tt]",
187 "a correlation plot of the Uij will be created, if any anisotropic",
188 "temperature factors are present in the pdb file.[PAR]",
189 "With option [TT]-dir[tt] the average MSF (3x3) matrix is diagonalized.",
190 "This shows the directions in which the atoms fluctuate the most and",
193 static bool bRes
=FALSE
,bAniso
=FALSE
,bdevX
=FALSE
,bFit
=TRUE
;
195 { "-res", FALSE
, etBOOL
, {&bRes
},
196 "Calculate averages for each residue" },
197 { "-aniso",FALSE
, etBOOL
, {&bAniso
},
198 "Compute anisotropic termperature factors" },
199 { "-fit", FALSE
, etBOOL
, {&bFit
},
200 "Do a least squares superposition before computing RMSF. Without this you must make sure that the reference structure and the trajectory match." }
202 int step
,nre
,natom
,natoms
,i
,g
,m
,teller
=0;
203 real t
,lambda
,*w_rls
,*w_rms
;
209 t_atoms
*pdbatoms
,*refatoms
;
214 int status
,npdbatoms
,res0
;
219 FILE *fp
; /* the graphics file */
220 const char *devfn
,*dirfn
;
228 real bfac
,pdb_bfac
,*Uaver
;
232 double *rmsf
,invcount
,totmass
;
239 char *leg
[2] = { "MD", "X-Ray" };
242 { efTRX
, "-f", NULL
, ffREAD
},
243 { efTPS
, NULL
, NULL
, ffREAD
},
244 { efNDX
, NULL
, NULL
, ffOPTRD
},
245 { efPDB
, "-q", NULL
, ffOPTRD
},
246 { efPDB
, "-oq", "bfac", ffOPTWR
},
247 { efPDB
, "-ox", "xaver", ffOPTWR
},
248 { efXVG
, "-o", "rmsf", ffWRITE
},
249 { efXVG
, "-od", "rmsdev", ffOPTWR
},
250 { efXVG
, "-oc", "correl", ffOPTWR
},
251 { efLOG
, "-dir", "rmsf", ffOPTWR
}
253 #define NFILE asize(fnm)
255 CopyRight(stderr
,argv
[0]);
256 parse_common_args(&argc
,argv
,PCA_CAN_TIME
| PCA_CAN_VIEW
| PCA_BE_NICE
,
257 NFILE
,fnm
,asize(pargs
),pargs
,asize(desc
),desc
,0,NULL
,
260 bReadPDB
= ftp2bSet(efPDB
,NFILE
,fnm
);
261 devfn
= opt2fn_null("-od",NFILE
,fnm
);
262 dirfn
= opt2fn_null("-dir",NFILE
,fnm
);
264 read_tps_conf(ftp2fn(efTPS
,NFILE
,fnm
),title
,&top
,&ePBC
,&xref
,NULL
,box
,TRUE
);
265 snew(w_rls
,top
.atoms
.nr
);
267 fprintf(stderr
,"Select group(s) for root mean square calculation\n");
268 get_index(&top
.atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),1,&isize
,&index
,&grpnames
);
271 for(i
=0; i
<isize
; i
++)
272 w_rls
[index
[i
]]=top
.atoms
.atom
[index
[i
]].m
;
274 /* Malloc the rmsf arrays */
277 for(i
=0; i
<isize
; i
++)
284 get_stx_coordnum(opt2fn("-q",NFILE
,fnm
),&npdbatoms
);
287 init_t_atoms(pdbatoms
,npdbatoms
,TRUE
);
288 init_t_atoms(refatoms
,npdbatoms
,TRUE
);
289 snew(pdbx
,npdbatoms
);
290 /* Read coordinates twice */
291 read_stx_conf(opt2fn("-q",NFILE
,fnm
),title
,pdbatoms
,pdbx
,NULL
,NULL
,pdbbox
);
292 read_stx_conf(opt2fn("-q",NFILE
,fnm
),title
,refatoms
,pdbx
,NULL
,NULL
,pdbbox
);
294 pdbatoms
= &top
.atoms
;
295 refatoms
= &top
.atoms
;
297 npdbatoms
= pdbatoms
->nr
;
298 snew(pdbatoms
->pdbinfo
,npdbatoms
);
299 copy_mat(box
,pdbbox
);
303 sub_xcm(xref
,isize
,index
,top
.atoms
.atom
,xcm
,FALSE
);
305 natom
= read_first_x(oenv
,&status
,ftp2fn(efTRX
,NFILE
,fnm
),&t
,&x
,box
);
307 /* Now read the trj again to compute fluctuations */
311 /* Remove periodic boundary */
312 rm_pbc(&(top
.idef
),ePBC
,natom
,box
,x
,x
);
314 /* Set center of mass to zero */
315 sub_xcm(x
,isize
,index
,top
.atoms
.atom
,xcm
,FALSE
);
317 /* Fit to reference structure */
318 do_fit(natom
,w_rls
,xref
,x
);
321 /* Calculate Anisotropic U Tensor */
322 for(i
=0; i
<isize
; i
++) {
324 for(d
=0; d
<DIM
; d
++) {
325 xav
[i
*DIM
+ d
] += x
[aid
][d
];
327 U
[i
][d
*DIM
+ m
] += x
[aid
][d
]*x
[aid
][m
];
332 /* Calculate RMS Deviation */
333 for(i
=0;(i
<isize
);i
++) {
335 for(d
=0;(d
<DIM
);d
++) {
336 rmsd_x
[i
][d
] += sqr(x
[aid
][d
]-xref
[aid
][d
]);
342 } while(read_next_x(oenv
,status
,&t
,natom
,x
,box
));
345 invcount
= 1.0/count
;
348 for(i
=0; i
<isize
; i
++) {
350 xav
[i
*DIM
+ d
] *= invcount
;
352 for(m
=0; m
<DIM
; m
++) {
353 U
[i
][d
*DIM
+ m
] = U
[i
][d
*DIM
+ m
]*invcount
354 - xav
[i
*DIM
+ d
]*xav
[i
*DIM
+ m
];
355 Uaver
[3*d
+m
] += top
.atoms
.atom
[index
[i
]].m
*U
[i
][d
*DIM
+ m
];
357 totmass
+= top
.atoms
.atom
[index
[i
]].m
;
359 for(d
=0; d
<DIM
*DIM
; d
++)
363 for(d
=0; d
<DIM
*DIM
; d
++) {
364 average_residues(NULL
,U
,d
,isize
,index
,w_rls
,&top
.atoms
);
369 for(i
=0; i
<isize
; i
++) {
371 pdbatoms
->pdbinfo
[aid
].bAnisotropic
= TRUE
;
372 pdbatoms
->pdbinfo
[aid
].uij
[U11
] = 1e6
*U
[i
][XX
*DIM
+ XX
];
373 pdbatoms
->pdbinfo
[aid
].uij
[U22
] = 1e6
*U
[i
][YY
*DIM
+ YY
];
374 pdbatoms
->pdbinfo
[aid
].uij
[U33
] = 1e6
*U
[i
][ZZ
*DIM
+ ZZ
];
375 pdbatoms
->pdbinfo
[aid
].uij
[U12
] = 1e6
*U
[i
][XX
*DIM
+ YY
];
376 pdbatoms
->pdbinfo
[aid
].uij
[U13
] = 1e6
*U
[i
][XX
*DIM
+ ZZ
];
377 pdbatoms
->pdbinfo
[aid
].uij
[U23
] = 1e6
*U
[i
][YY
*DIM
+ ZZ
];
385 for(i
=0; i
<isize
; i
++)
386 rmsf
[i
] = U
[i
][XX
*DIM
+ XX
] + U
[i
][YY
*DIM
+ YY
] + U
[i
][ZZ
*DIM
+ ZZ
];
389 fprintf(stdout
,"\n");
390 print_dir(stdout
,Uaver
);
391 fp
= ffopen(dirfn
,"w");
396 for(i
=0; i
<isize
; i
++)
400 /* Write RMSF output */
402 bfac
= 8.0*M_PI
*M_PI
/3.0*100;
403 fp
= xvgropen(ftp2fn(efXVG
,NFILE
,fnm
),"B-Factors",
404 label
,"(A\\b\\S\\So\\N\\S2\\N)",oenv
);
405 xvgr_legend(fp
,2,leg
,oenv
);
406 for(i
=0;(i
<isize
);i
++) {
407 if (!bRes
|| i
+1==isize
||
408 top
.atoms
.atom
[index
[i
]].resind
!=top
.atoms
.atom
[index
[i
+1]].resind
) {
409 resind
= top
.atoms
.atom
[index
[i
]].resind
;
410 pdb_bfac
= find_pdb_bfac(pdbatoms
,&top
.atoms
.resinfo
[resind
],
411 *(top
.atoms
.atomname
[index
[i
]]));
413 fprintf(fp
,"%5d %10.5f %10.5f\n",
414 bRes
? top
.atoms
.resinfo
[top
.atoms
.atom
[index
[i
]].resind
].nr
: i
+1,rmsf
[i
]*bfac
,
420 fp
= xvgropen(ftp2fn(efXVG
,NFILE
,fnm
),"RMS fluctuation",label
,"(nm)",oenv
);
421 for(i
=0; i
<isize
; i
++)
422 if (!bRes
|| i
+1==isize
||
423 top
.atoms
.atom
[index
[i
]].resind
!=top
.atoms
.atom
[index
[i
+1]].resind
)
424 fprintf(fp
,"%5d %8.4f\n",
425 bRes
? top
.atoms
.resinfo
[top
.atoms
.atom
[index
[i
]].resind
].nr
: i
+1,sqrt(rmsf
[i
]));
429 for(i
=0; i
<isize
; i
++)
430 pdbatoms
->pdbinfo
[index
[i
]].bfac
= 800*M_PI
*M_PI
/3.0*rmsf
[i
];
433 for(i
=0; i
<isize
; i
++)
434 rmsf
[i
] = (rmsd_x
[i
][XX
]+rmsd_x
[i
][YY
]+rmsd_x
[i
][ZZ
])/count
;
436 average_residues(rmsf
,NULL
,0,isize
,index
,w_rls
,&top
.atoms
);
437 /* Write RMSD output */
438 fp
= xvgropen(devfn
,"RMS Deviation",label
,"(nm)",oenv
);
439 for(i
=0; i
<isize
; i
++)
440 if (!bRes
|| i
+1==isize
||
441 top
.atoms
.atom
[index
[i
]].resind
!=top
.atoms
.atom
[index
[i
+1]].resind
)
442 fprintf(fp
,"%5d %8.4f\n",
443 bRes
? top
.atoms
.resinfo
[top
.atoms
.atom
[index
[i
]].resind
].nr
: i
+1,sqrt(rmsf
[i
]));
447 if (opt2bSet("-oq",NFILE
,fnm
)) {
448 /* Write a pdb file with B-factors and optionally anisou records */
449 for(i
=0; i
<isize
; i
++)
450 rvec_inc(xref
[index
[i
]],xcm
);
451 write_sto_conf_indexed(opt2fn("-oq",NFILE
,fnm
),title
,pdbatoms
,pdbx
,
452 NULL
,ePBC
,pdbbox
,isize
,index
);
454 if (opt2bSet("-ox",NFILE
,fnm
)) {
455 /* Misuse xref as a temporary array */
456 for(i
=0; i
<isize
; i
++)
458 xref
[index
[i
]][d
] = xcm
[d
] + xav
[i
*DIM
+ d
];
459 /* Write a pdb file with B-factors and optionally anisou records */
460 write_sto_conf_indexed(opt2fn("-ox",NFILE
,fnm
),title
,pdbatoms
,xref
,NULL
,
461 ePBC
,pdbbox
,isize
,index
);
464 correlate_aniso(opt2fn("-oc",NFILE
,fnm
),refatoms
,pdbatoms
,oenv
);
465 do_view(oenv
,opt2fn("-oc",NFILE
,fnm
),"-nxy");
467 do_view(oenv
,opt2fn("-o",NFILE
,fnm
),"-nxy");
469 do_view(oenv
,opt2fn("-od",NFILE
,fnm
),"-nxy");